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PREFACE 

The  results  reported  herein  were  obtained  from  the 
Arnold  Engineering  Development  Center  by  the  University  of 
Cincinnati,  Department  of  Aerospace  Engineering  under  Con- 
tract F40600-74-C-0011.  The  authors  of  this  report  were 
B.  N.  Srivastava,  M.  J.  Werle,  and  R.  T.  Davis.  The  Air 
Force  Project  Engineer  for  this  contract  was  E.  R.  Thompson, 
AEDC/DYR.  The  Program  Element  Number  was  65807F. 

This  report  covers  the  work  during  the  period  March 
1976  to  October  1976.  The  reproducibles  used  in  the  repro- 
duction of  this  report  were  supplied  by  the  authors. 
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1.  INTRODUCTION 

Recent  developments  in  aerodynamics  and  space  flight 
have  caused  increasingly  focused  attention  on  the  problem 
of  theoretically  predicting  the  blunt  body  flow  field. 

Three  current  numerical  approaches  for  treating  this  problem 
include  solution  of  either  the  full  Navier-Stokes  equations, 
the  second-order  boundary-layer  equations,  or  the  viscous 
shock  layer  equations.  Use  of  the  full  Navier-Stokes 
equations  [1]  has  heen  quite  successful  in  providing 
solutions  for  stagnation  regions,  but  generally  have  been 
applied  for  only  one  nose  radius  downstream.  This  is 
because  the  elliptic  nature  of  these  equations  increases 
the  complexity  of  the  solution  procedure  and  restricts 
the  application  of  these  methods  in  the  downstream 
direction.  While  there  are  several  computational  diffi- 
culties associated  with  the  second-order  boundary  layer 
approach  [2],  many  of  the  difficulties  associated  with 
computing  viscous  hypersonic  flows  over  blunt  bodies  can 
be  overcome  through  use  of  the  viscous  shock  layer 
equations.  In  this  approach  the  entire  flow  field,  from 
the  body  to  the  shock,  is  treated  in  a unified  manner. 

This  is  the  approach  taken  here,  wherein  the  basic  method 
of  Davis  [2]  is  applied  to  nonanalytic  bodies  such  as 
spherically  blunted  cones. 
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Blunted  cones  merit  further  study,  both  because  some 
blunting  of  the  tip  of  a sharp  cone  due  to  extreme 
heating  appears  unavoidable  and  because  blunting  has  favor 
able  effects  on  transition.  Previous  research  [3,  4]  has 
shown  that  a small  amount  of  bluntness  to  a sharp  cone  is 
conducive  to  a delay  in  boundary  layer  transition.  Unfor- 
tunately only  limited  work  has  been  done  in  the  past  for 

■ 

such  nonanalytic  bodies,  in  part,  due  to  the  difficulties 
associated  with  the  discontinuity  occurring  in  the 
surface  curvature  at  the  sphere/cone  juncture  point  when- 
ever a surface  coordinate  system  is  used.  Miner  and 
Lewis  [5]  and  Kang  and  Dunn  [6]  have  been  able  to  obtain 
some  numerical  solutions  for  such  nose  cone  problems. 

The  approach  of  Miner  and  Lewis  [5]  was  to  smoothen  the 

> 

effect  of  the  curvature  discontinuity  at  the  sphere/cone 
tangency  point  by  constructing  an  artificially  continuous 
distribution  of  curvature.  Kang  and  Dunn  [6]  approached 
the  same  problem  by  application  of  a Karman-Pohlhausen 
integral  method  to  the  Navier-Stokes  equations  under  a 
thin  shock  layer  assumption.  They  treated  the  nose  region 
as  well  as  the  afterbody  conical  section  of  the  spheri- 
cally blunted  cones  in  a similar  manner.  While  these 
solutions  yield  acceptable  numerical  results , it  appears 
that  a more  consistant  formulation,  which  accounted  for 
the  flow  behavior  at  the  sphere  cone  tangency  point,  could 
improve  these  solutions  and  enhance  their  reliability. 
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The  basic  difficulty,  however,  is  that  only  a very 
limited  amount  of  information  is  available  for  the  flow 
behavior  at  such  a curvature  discontinuity  - information 
that  could  be  utilized  in  a computational  procedure  to 
achieve  more  accurate  solutions.  It  is  known  from 
inviscid  considerations  that  a streamline  curvature 
discontinuity  produces  a discontinuity  in  the  flow 
gradients  along  that  streamline.  Viscous  effects  in  such 
problems  are  expected  to  remove  such  surface  gradient 
jumps  in  a manner  shown  by  Messiter  and  Hu  [7]  for  two 
dimensional  flows.  Their  analysis  for  high  Reynolds 
number  flows  shows  that  the  jump  in  pressure  gradient 
predicted  by  the  inviscid  flow  theory  is  smoothed  by 
viscous  effects  through  a triple  deck  structure  near  the 
juncture  region.  However  the  region  within  which  this 
smoothing  takes  place  is  found  to  be  small  for  realistic 
finite  Reynolds  numbers  and  it  would  be  anticipated  that 
a viscous  numerical  calculation  with  a finite  mesh  size 
might  not  capture  this  physical  phenomenon  and  would 
therefore  still  predict  a gradient  discontinuity.  Thus, 
one  is  faced  with  the  prospect  of  generating  finite 
difference  approximations  for  regions  in  which  the  flow 
gradients  are  virtually  discontinuous. 

An  added  complication  arises  because  the  viscous 
shock  layer  equations,  when  written  in  a surface  coordinate 
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system,  contain  explicit  dependence  on  the  surface  curvature. 
The  surface  curvature  undergoes  a discontinuous  change  from 
a value  of  one  on  the  sphere  to  a value  of  zero  on  the  cone 
portion  resulting  in  the  appearance  of  gradient  discon- 
tinuities in  all  flow  variables  at  all  points  normal  to 
the  body  surface  across  the  shock  layer  thickness  at  the 
juncture  point.  These  discontinuities  are  purely  an 
artifice  and  occur  only  because  of  the  choice  of  surface  1 
coordinates.  None-the-less  they  must  be  accounted  for  in 
any  numerical  approximating  procedure  if  reliable  results 
are  to  be  achieved. 

An  analytical  assessment  of  the  flow  behavior  at  the 

f 

sphere/cone  juncture  point  and  subsequently  a proper  numerical 
formulation  of  this  problem  is  the  purpose  of  the  present 
study.  It  is  believed  that  when  a surface  coordinate 
frame  of  reference  is  used  a finite  difference  formulation 
of  the  governing  equations  must  be  such  that  the  longitudinal 
derivatives  be  carefully  evaluated  at  the  sphere/cone 
juncture  point  in  order  to  eliminate  large  numerical 
truncation  errors.  This  can  be  done  by  ensuring  that  the 
finite  difference  form  of  the  longitudinal  derivatives 
avoid  any  differencing  across  the  sphere/cone  juncture 
discontinuity.  This  technique  is  demonstrated  in  the 
present  approach  where  the  numerical  scheme  utilizes  a 
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time  dependent  relaxation  technique  for  the  shock  wave 
shape.  A model  problem  analogous  to  the  sphere/cone 
juncture  problem  is  first  formulated  and  finite  difference 
schemes  developed . and  demonstrated  for  this  case.  It  is 
shown  that  the  present  method  accurately  captures  the 
anticipated  discontinuous  behavior  of  the  flow  derivatives 
at  the  model  juncture  point.  This  concept  is  then 
extended  for  the  solution  of  the  viscous  shock  layer 
equations  for  hypersonic  flow  past  spherically  blunted 
cones  with  half  cone  angles  varying  from  30°  to  0°  under 
various  free  stream  conditions.  The  results  indicate 
good  comparisons  with  inviscid  solutions  and  experimental 
data. 
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II.  GOVERNING  EQUATIONS 

The  viscous  shock  layer  concept  has  been  presented  in 
detail  by  Davis  [2]  and  therefore  is  only  summarized  here. 
The  compressible  Navier-Stokes  equations  are  written  in  a 
boundary  layer  like  coordinate  system  (see  Fig.  1)  and 
nondimensionalized  by  variables  which  are  of  order  one  in 
the  region  near  the  body  surface  (boundary  layer)  for 
large  Reynolds  numbers.  The  same  set  of  equations ’ are 
then  written  in  variables  which  are  of  order  one  in  the 
essentially  inviscid  region  outside  the  boundary  layer. 

In  the  final  set  of  equations  terms  are  retained  that  are 
second-order  in  the  inverse  square  root  of  a Reynolds 
number.  A comparison  of  the  two  sets  of  equations  is 
then  made  and  one  set  of  equations  is  found  from  them 
which  is  valid  to  second  order  in  both  the  (inviscid) 
outer  and  inner  (viscous)  regions.  A solution  to  this 
set  of  equations  is  thus  uniformly  valid  to  second  order 
in  the  entire  shock  layer  for  arbitrary  Reynolds  number. 

i 

The  resulting  equations  (and  notation)  are  the  same  as 
those  presented  by  Davis  [2]  and  are  given  as: 

Continuity: 

[ (r+ncos$)  ^pu]  + [ (l+«n)  (r+ncos<j>)  ^pv]  „ = 0 (la) 

£>  7i 
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Longitudinal  Momentum: 

p{u  u /(1+Kn)  + v u + Kuv/{l+icn)}  + p /(1+Kn) 

S I*  s 

= [e2/  (l+icn) 2 (r+ncos<|>)  ■*]  [ (l+xn)  2 (r+ncos<|>)  ^ i]n 


(lb) 

where, 

t = y tun  ” ieu/(l+icn)]  ' (lc) 

Normal  Momentum: 

2 

p{u  v_/(l+icn)  + v v - ku  /(1+Kn)}  + p = 0 (Id) 

s n n 

which  with  the  thin  shock  layer  approximation  becomes, 

-p  ku2/ (1+Kn)  + pn  = 0 (le) 

Energy  Equation: 

p {uTg/ (1+Kn)  + vTnJ-  u ps/(l+Kn)  " v Pn  “ e2t2/m 
+ [E2/(1+Kn)  (r+ncos<Ji)  -*]  [(1+Kn)  (r+ncosij>)  ^qln  (If) 


where 


q = u T n/o  (lg) 

Equation  of  State: 

p = (y-DpT/y  (lh) 

Viscosity  Law: 

p = T3/2(l+c’)/(T+c')  (li) 
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where 

» * 2 * 
c = c /M„  T„  (y-D 

* 

and  c is  taken  to  be  198. 6 °R  for  air. 

The  boundary  conditions  employed  here  are  the  no 
slip  surface  conditions, 

u(s,o)  = v(s,o)  = 0 (2a) 

and 

T(s,o)  = Tw  (2b) 

while  at  the  shock  location  the  oblique  shock  relations 
are  used  to  relate  the  flow  variables  just  aft  of  the 
shock  to  the  free  stream  conditions  through  the  local 
shock  slope.  These  relations  are  given  as  ' 

ush  = ush  sin(a+&>  + vsh  cos (a+8 ) (3a) 

vsh  = - ugh  cos ( a+3 ) + vgh  sin  (a+fi)  (3b) 

where  ugh  and  vgh  are  velocity  components  at  the  shock 
given  as 

ugh  = cosa  (3c) 

vsh  = - sin a/p sh  (3d) 

and 

psh  = ^psh/(lf-1)Tsh  <3e> 
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with 

psh  = [2/  (Y+l)  1 sin2a  - CY-1)/Y  (Y+l)Mf  (3f) 

Tsh  = ^ush  “ cosa)2/2  + { [4y/ (y+l) 2] sin2a+ [2/ (y-1) 

-4(y-l)/(y+l)2)/M2  - 4/(Y+l)2Mf  sin2a}/2  (3g) 

The  shock  angle,  a,  is  related  to  the  shock  thickness, 
ng,  through  the  relation 

cin 

/ (1+Kng)  = tan (a-<J>)  (4a) 

The  value  of  n„  itself  can  be  written  from  mass  conservation 
s 

considerations  as 

{r+n  cos<f>)1+^  = 2?  I pu  (r+ncosi)>)  ^ dn  (4b) 

s o 

For  reasons  explained  in  Reference  [2] , the  above 
equations  will  be  normalized  according  to  the  following 
scheme: 


" = n/nsh 

i = s 

u = u/ush 

v - v/vsh 

t = T/Tsh 

P - P/Psh 

* = t‘/<’sh 

S ■ v/vsh 

(5a-h) 

The  differential  relations  needed  to  transform  equations 
(la)  through  (li)  are  given  by 
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3/3 s = 3/3E  ~ nCnsh/ngh)  3/3n 


3/3n  = (l/ngh)  3/3n 

32/3n2  = d/n2h)  32/3n2  (6) 


where , 


n 


sh 


(ansh/dC) 


The  s-momentum  and  energy  equations  (lb)  and  (If)  written 
in  the  transformed  E,n  plane  can  be  conveniently  put  in 
a standardized  form  for  a parabolic  equation  as, 


2 2 

3 w/3n  + c^Ow/Sn)  + ct2w  + a3  + a4(3w/3E)  =0  (7) 

where  w represents  u for  the  s-momentum  equation,  and 
t for  the  energy  equation.  For  the  momentum  equation  the 
coefficients  ct^  ->  are: 


pshushnsh 

al 1 

e ysh 


n 

1+Kn 


sh 

sh 


n 


pun 

v 


pshvshnsh 

2 

e psh 


pv 

u 


, - . Knsh 

+ V1*  + T+Z?r 


cos<j>n 


sh 


shn  Y+nshncos^> 


(8a) 
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“2  = 


pshushnsh  nsh  £u  _ pshvshnBh  ‘ Knsh  - 

1+Kn  un  - 2 l+ien  , n p 

sh  y e y _i_  sh 


2 

e M 


sh 


sh 


n.  kii 

K iraw?  V"  * Ctt_ 


sh  . °°s*nsh 


-)  x(. 


<n 


sh 


■l.+lcnshn  r+ns^ncos(j)  1+icn 


(8b) 


°3  = 


pshnsh 

“2 

e ysh 


n 


sh 


n 


sh  - 


sh 


l/u,h(PE~«P„^P) 


sh 


sh 
(8c) 


“4  = ' (',shushnsh/£;2’,sh)  <nsh/(1+ltnsh,,)  5s 


(8d) 


for  the  energy  equation  the  coefficients  are 


_ pshushnsha  nsh  pun  _ pshvshnsha  pv 

“1 2 l+raTT  2 x - 


e y 


sh 


■kii  , n 

sh  y e y 


sh 


- - Knsh  cos^nsh 

+ y /y  + — + — — 

pn  ^+Knshn  r+ng^ncosifr 


(9a) 


a2  ” ~(pshushTsh/e  ^sh^h*  (nsh/(1+Knshn))  x -U  (9b) 

y 


a3  2 


pshushnsh° 


e ^sl^sh 


- nshU  - nsh  - psh  - 

V-  ^ <p? + 


, Sh  -r  . w‘shu  Knsh  -,2 

+ vp"1+  *S"  <u"  u) 


(9c) 
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“4  “ ■tapshushnsh/£2“Sh)  (nsh/(1+Knshn)  )^ 


<9d) 


The  remaining  differential  equations  are  first  order  and 
are  the  continuity  equation: 


[nsh(r+nshncos<|,)  pshush  pG]5  + C(r+nshncos<fr)x 


{(1+<nshn)PshVsh  pv  - nsh  pshush  pun}] ^ = 0 


(10) 


and  the  n-momentum  equation : 

1 

<1+K"sh,|)  <V?  ' nsh/nah  nV,>  + ^ V>+  “Z  V,> 


u 


1+Knsh"  vsh 


sh  --2  , 
pu  + 


sh 


p t-U  , v ,n  , *■  n 
Ksh  sh  sh  sh 


P = 0 


(11) 


where  with  the  thin  layer  approximation  this  equation 
becomes, 

P„  ■ («/<1+l'nsh"1 1 <pshushnsh/psh)  p"2 


(12) 


This  leaves  the  equation  of  state. 


p = pt 


(13) 


and  the  viscosity  law, 


V = I(Tgh  + c )/(Tght  + c )]  t3/2 


(14) 
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At  the  shock  location  all  variables  are  unity, 

u = v = t = p=  p = v = l at  n = 1 (15) 

An  equation  of  mass  conservation  can  be  obtained  from 
equation  (10)  by  integrating  from  n = 0 to  n = 1 .-while 
holding  £ constant.  This  results  in 

af-  <r+nshoos*)|nsh‘,Shush  ■(1+Knsh,'IShvsh1  (16a) 

where 

1 

m = I nsh (r+nshnCOS*J  pshush  dn  * 

o 

\ 

is  proportional  to  the  rate  of  mass  flux  between  the  body  n 
and  shock  at  a given  position  on  the  body  surface. 

Equations  (7),  (10)-(14)  and  (16),  constitute  the 
complete  set  of  governing  equations  for  the  unknowns  u,  v, 
t,  p,  y,  p and  n^.  These  equations  are  solved  along 
with  the  surface  boundary  conditions  given  by  equations 
(2a)  and  (2b)  and  the  shock  conditions  given  by  equation 
(15) . The  mass  conservation  equation  (16a)  and  (16b)  is 
used  to  determine  the  shock  stand  off  distance  n^.  The 
general  procedure  is  to  evaluate  the  rate  of  mass  flux 
between  the  body  and  shock  at  a given  position  on  the 
body  surface  from  equation  (16b)  using  the  values  of  the 
physical  quantities  previously  calculated  and  then  evaluate 
n^  from  equation  (16a)  . 
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III.  JUNCTURE  REGION 
1.  Inviscid  Flow  Analysis 

The  aim  of  the  present  section  is  to  analyze  the  nature 
of  the  flow  at  a sphere/cone  juncture  point  in  the  light  of 
various  forms  of  the  gas  dynamic  equations.  It  is  of 
interest  to  analyze  the  viscous  as  well  as  inviscid  gas 
dynamic  equations  in  order  to  understand  the  physical 
behavior  at  a point  in  the  flow  field  where  a discontinuity 
in  curvature  is  encountered.  In  order  to  be  fully  con- 
sistent in  the  present  analysis,  it  is  desirable  to  first 
consider  the  flow  behavior  from  an  inviscid  standpoint  and 
then  subsequently  include  viscous  effects. 

From  an  inviscid  standpoint,  it  is  known  from  Euler* s‘ 
equations  that  a streamline  curvature  discontinuity 
produces  a discontinuity  in  the  flow  gradients  only  along 
that  streamline  [8].  For  locally  supersonic  flows,  the 
discontinuity  at  a juncture  point  (e.g.  sphere/cone  tan- 
gency  point)  in  the  flow  gradients  will  be  propagated  along  the 
characteristic  lines  inclined  at  the  local  Mach  angle  to 
the  flow  direction.  Eventually,  therefore,  such  discon- 
tinuities in  the  flow  gradients  would  propagate  everywhere 
within  the  flow  field  downstream  of  the  juncture  discontinuity 
due  to  successive  reflections  of  these  characteristic 
lines  from  the  shock  and  body  surface. 
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Note,  however,  that  at  the  juncture  point  the  flow  gra- 
dients will  be  continuous  across  the  shock  layer  except  at  the 
body  surface  where  gradient  discontinuities  will  be 
present..  Note  must  also  be  made  here  that  the  flow 
variables  themselves  are  found  to  be  continuous  at  the 
juncture  point  across  the  shock  layer.  However,  an 
added  complication  arises  due  to  the  explicit  appearance 
of  the  surface  curvature  in  these  equations  when  written 
in  a surface  coordinate  system.  Since  the  surface 
curvature  itself  is  discontinuous  at  the  sphere/cone 
juncture  point,  it  is  necessary  to  assess  the  influence 
of  the  discontinuity  on  the  flow  properties  and  their 
derivatives.  Intuitively  it  is  obvious  that  a mere 
coordinate  transformation  would  not  affect  the  physical 
behavior  so  that  the  flow  variables  themselves  are 
continuous  at  the  juncture  point  all  across  the  shock 
layer.  This  can  also  be  straightforwardly  demonstrated 
through  consideration  of  the  integral  form  of  the  con- 
servation laws  in  the  surface  coordinate  system  as  shown 
in  Appendix  (A) . 

Note  that  the  set  of  equations  (A 7 , 14,  20)’  in  Appendix 
(A)  provide  for  either  the  trivial  case  of  p^  = P2  and 
U1  = u2'  or  a shock  like  discontinuity.  Since  there  is  no 
physical  event  that  could  cause  a shock  at  the  sphere/cone 
juncture  point,  it  is  fairly  evident  that  the  trivial 


21 


AE  DC-TR-77-20 


solution  is  expected  for  this  case  indicating  that  the 
flow  variables  are  continuous  at  the  sphere/cone  juncture 
point  in  the  surface  coordinate  system.  However  the 
same  conclusion  does  not  apply  to  the  flow  derivatives 
with  respect  to  the  surface  distance.  This  is  evident 
from  the  differential  form  of  the  two-dimensional  inviscid 
gas  dynamic  equations  as  recovered  from  the  integral 
equations.  These  are  given  as. 

Continuity 

(pu)g  + (l+<n)  ( p v) n + kPv  = 0 (17) 

s -Momentum 

p_  + puu  + (1+icn)  pvu„  + k p'uv  = 0 (18) 

o S lx 

n -Momentum 

2 

puv_  /(l+tcn)  + pvvn  - k p u /(1+icn)  + p„  = 0 (19) 

s n n 

Energy  Equation 

puT  + pvT  (l+<n)  - up  - vp  (1+icn)  = 0 (20) 

s v n s n 

Equation  of  State 

p = (III)  PT  (21) 

2 

Using  the  equation  of  state  and  defining’  a = lyp/p),  the 
energy  equation  can  be  rewritten  as, 

2 2 

up  + v(l+Kn)p_  - ua  p - v(l+icn)a  = 0 (22) 

5 ns  n 
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Equations  (17)  through  (22)  can  now  be  used  on  the  two 
sides  of  the  sphere/cone  juncture  point,  noting  that  the 
surface  curvature,  k,  takes  a value  of  1.0  on  the  spherical 
part  and  a value  of  0 on  the  conical  part.  Using  subscript 
1 and  2 for  the  sphere  and  cone  portions  respectively 
and  noting  that  the  flow  variables  are  continuous  across 
the  juncture  point,  the  inviscid  equations  provide  the 
following  jump  conditions  at  the  juncture  point, 

(23a) 
(23b) 

(23c) 
(23d) 

In  addition  to  the  flow  variables  themselves,  consi- 
deration must  also  be  given  to  the  bow  shock  shape.  Note 
that  for  problems  of  this  type  the  shock  shape  itself  would 
be  expected  to  be  smooth  through  the  sphere  cone  juncture 
point  with  the  shock  slope  at  any  location  given  in  terms 
of  the  axial  coordinate  system  (Figure  1)  as 

tana  = (24) 


ps  = (1+n>  Ps 
S1  S2 


p = (1+n)  p 
S1  S2 


u = (1+n)  u — v 


s 


v = (1+n)  v + u 


s. 


23 


AEDC-T  R-77-20 


Since  the  shock  shape,  R,  itself  is  a smooth  function  through 
the  sphere/cone  tangency  point,  the  shock  slope,  a, 
would  be  smooth  at  this  point.  However,  the  first  derivative 
of  the  shock  distance,  R,  with  respect  to  the  surface 
coordinate  system  is  obtained  from  the  geometrical  relation; 


dR 

ds 


Cl+Kn  ) 
s 


sing 

cos  ta-$) 


(25) 


Since  the  shock  angle,  a,  and  the  body  angle,  <J>,  are 
continuous  functions  of  the  surface  distance,  this  relation 
yields  a jump  condition  for  dR/ds  at  the  sphere/cone 
tangency  point  as, 


sphere 


■ (1+ns><af> 


cone 


(26) 


Similar  discontinuous  behavior  can  be  shown  to  appear  in 
derivatives  such  as  dn  /ds  and  dx  /ds  as  shown  by  the 

S S 

following  expressions, 


(1+icn  ) tan(a-ij)) 

5 


(27) 


dx 

£ 

ds 


= (l+lcn  ) .co.^ 

s cos  (g-<f>) 


(28) 


The  corresponding  jump  conditions 


written  as, 
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dns 

as  sphere 


dn 


tl+ns>  izr* 


cone 


(29) 


dx. 


dx. 


(di£)  v = (1+ns)  (di^) 
u“  sphere  s as  cone 

It  is  of  interest  to  note  that  the  expression  (24) , 
rewritten  using  the  surface-coordinate  system  as, 


(30) 
can  be 


dx 

ds 


s 


tana 


dR 

ds. 


indicating  that  the  discontinuity  associated  with  dxs/ds 
and  dR/ds  at  the  juncture  point  are  of  the  nature  such 
that  the  shock  slope,  a,  itself  is  continuous  at  this 
point. 

It  must  be  pointed  out  here  that  similar  jump  conditions 
can  also  be  established  for  the  second  derivatives  of  the 
flow  quantities  mentioned  above,  whenever  required.  These 
derivatives  are  undefined  at  the  juncture  point  in  this 
coordinate  system,  however  finite  values  exist  for  these 
quantitites  immediately  ahead  and  behind  the  juncture  point. 

A typical  case  where  higher  derivatives  are  needed  is  at 
the  shock  location.  The  derivatives  of  the  flow  properties 
behind  the  shock  are  shown  in  Appendix  (B) . Note  that 
these  derivatives  (Bll,  B 15-17)  undergo  discontinuous 
changes  at  the  sphere/cone  juncture  point  and  the 
magnitude  of  these  discontinuities  are  related  to  the 
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surface  curvature,  k,  and  the  discontinuity  associated 

2 2 

with  the. second  derivative  of  the  shock  shape,  d R/ds  . 
However  it  is  important  to  emphasize,  here,  that  this 
situation  is  entirely  due  to  the  use  of  the  surface 
coordinate  system. 

It  is,  thus,  found  that  within  the  framework  of 
Euler's  equations  the  flow  properties  are  continuous 
whereas  the  flow  derivatives  with  respect  to  surface 
distance  are  discontinuous  across  the  layer  at  the 
sphere/cone  juncture  point.  The  slope  of  the  shock 
relative  to  the  surface  distance  is  also  found  to  be 
discontinuous  at  this  point. 

2.  Viscous  Flow  Analysis 

(a)  Classical  Boundary  Layer 

From  a viscous  standpoint  it  seems  more  rational  to 
first  address  the  question  of  validity  of  the  various  forms 
of  the  viscous  gas  dynamic  equations  as  applied  to  the 
spherically  blunted  cones  with  a discontinuous  surface 
curvature  at  the  juncture  point.  The  boundary  layer 
version  of  the  Navier-Stokes  equations  cannot  hold  at  the 
sphere/cone  juncture  point  becuase  the  longitudinal 
derivative  of  the  surface  curvature,  3k/3s  >>  1 [9,  10] 

and  thus  the  gradient  of  the  corresponding  inviscid  ■ 
pressure  is  discontinuous  there.  Any  such  discontinuity 
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would  seem  to  be  in  violation  of  the  boundary  layer  scaling 
laws  wherein  longitudinal  derivatives  are  assumed  to  be 
much  smaller  than  normal  derivatives.  This  becomes  more 
apparent  when  one  considers  inclusion  of  higher  order 
terms  in  the  boundary  layer  equations.  The  second  order 
correction  (in  Re”^2)  due  to  longitudinal  curvature  is 
driven  by  the  rate  of  change  of  curvature  - thus  causing 
this  higher  order  effect  to  rise  up  to  first  order  level 
near  a sphere/cone  juncture  point.  It  is>  therefore, 
apparent  that  a new  local  solution  needs  to  be  developed 
near  the  juncture  point  in  order  to  accommodate  this 
anomoly.  Such  an  analysis  has  been  performed  by  Messiter 
and  Hu  [7]  for  two-dimensional  flows. 

(b)  Triple  Deck  Analysis 

A study  of  the  juncture  region  has  been  completed  by 
Messiter  and  Hu  [7]  for  two-dimensional  flow  problems. 

They  point  out  that,  unlike  the  classical  boundary  layer 
case,  an  interaction  with  the  external  flow  must  be 
taken  into  account,  this  occurring  through  a small  pressure 
change  acting  over  a suitably  small  distance  along  the 
boundary  layer.  The  details  of  the  resulting  local 
pressure  distribution  cannot  be  specified  in  advance, 
but  must  be  found  by  studying  changes  in  the  boundary 
layer  coupled  with  small  perturbations  on  the  external 
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flow.  Messiter's  analysis  shows  that  the  discontinuity 
in  the  pressure  gradient  predicted  by  the  inviscid  flow 
theory  can  be  removed  by  using  a triple  deck  formulation 
and  continuous  expressions  for  the  pressure  gradient  can 
be  obtained  which  are  presumed  to  be  correct  asymptotic 
representations  as  the  viscosity  coefficient  approaches 
zero.  This  is  achieved  by  noting  that,  locally, the  most 
important  changes  in  the  profile  shape  occur  in  a thin 
sublayer  [11,  12]  close  to  the  wall  where  the  changes 
in  the  viscous,  pressure, and  inertia  forces  are  all  of  the 
same  order  as  the  characteristic  Reynolds  number  tends 
to  infinity.  The  remainder  of  the  boundary  layer 
experiences  primarily  a displacement  effect  because  of 
the  small  acceleration  of  the  fluid  in  the  sublayer,  and 
the  resulting  small  decrease  in  the  flow  deflection  angle 
is  nearly  constant  across  most  of  the  boundary  layer. 

The  interaction  of  the  boundary  layer  with  the  external 
flow  occurs  in  a streamwise  distance,  X = 0 (Re3//8)  ,and 
the  sublayer  thickness  is  given  by  Y = OtRe^^8),  while 
the  pressure  change  is  found  to  be  of  0(Re3//8).  The 
present  sphere/cone  problem  is  more  complex  due  to  the 
axisymmetric  nature  of  tfie  Body  and  the  fact  that  the 
approaching  boundary  layer  at  the  juncture  point  is  not 
that  due  to  a flat  plate  as  it  was  in  Messiter  and  Hu's 
analysis.  However,  an  approximate  calculation  can  be 
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performed  using  their  analysis  to  determine  the  scale 
within  which  the  viscous  smoothing  takes  place  at  the 
sphere/cone  juncture  point.  This  can  be  done  by  using 
the  local  Reynolds  number  at  the  sphere/cone  juncture 
point.  Figure  C2)  shows  the  results  of  such  an  analysis 
where  the  surface  pressure  is  shown  against  the  distance 
in  physical  coordinates.  The  asymptotic  smoothing  of 
the  inviscid  pressure  gradient  at  the  juncture  point  is 
seen  to  be  achieved  for  this  case  in  a very  small  physical 
distance  upstream  and,  comparatively  larger, but  yet 
small  distance  downstream. 

(c)  Viscous  Shock  Layer 

These  results  imply  two  important  points.  First,  near 
the  sphere/cone  juncture  region  the  correct  asymptotic 
solution  can  be  obtained  provided  the  viscous  set  of  gas 
dynamic  equations  are  such  that  they  retain  the  boundary 
layer  and  inertia  terms  in  the  viscous  region  and  allow 
for  displacement  interaction  with  the  inviscid  flow. 

One  way  to  ensure  this  criteria  is  to  use  the  full 
Navier-Stokes  equations.  However,  the  full  viscous  shock 
layer  equations  also  seem  to  be  sufficient  since  they 
contain  all  the  viscous  terms  of  the  triple  deck  model 
plus  the  inertia  terms  that  take  into  account  interaction 
effects  in  the  inviscid  flow. 
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The  second  important  point  is  that  the  interaction 
effects  will  be  significant  in  only  a very  small  region 
of  the  physical  flow  and  will  be  difficult  to  detect 
for  high  Reynolds  number  cases. 

However,  note  that  the  present  choice  of  the 
coordinate  system  would  introduce  discontinuities  in  the 
longitudinal  flow  gradients  at  the  sphere/cone  juncture 
point,  as  observed  for  the  inviscid  flow.  It  can  be  shown 
also  that  when  viscous  effects,  as  included  in  the  full 
shock  layer  equations,  are  accounted  for,  the  flow  variables 
themselves  are  continuous  through  the  sphere/cone  juncture 
point  (see  Appendix  A) . This  can  be  done  by  considering 
the  integral  form  of  the  viscous  equations  and  by  evaluating 
them  for  an  infinitesimally  small  element  in  the  surface 
coordinate  system  (Appendix  A) . Note  that  once  again, 
the  set  of  equations  (A 7 , 14 , 20)  give  either  a trivial 
solution  yielding  p^  = P2  and  u^  = Uj  or  a shock  like 
jump  discontinuity.  It  is  observed  that  since  there  is 
no  physical  event  that  could  cause  a shock  at  the  sphere/ 
cone  juncture  point,  the  trivial  solution  is  the  only 
possibility  indicating  that  the  flow  variables  are 
continuous  at  the  juncture  point  for  the  full  shock  layer 
equations  in  the  surface  coordinate  system.  However  the 
same  conclusion  does  not  apply  to  the  flow  derivatives 
with  respect  to  the  surface  distance.  This  is  evident  from 
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the  differential  form  of  the  full  shock  layer  equations. 
These  full  shock  layer  equations  (la-h)  can  be  used 
to  determine  the  jump  conditions  on  the  two  sides  of 
the  sphere/cone  juncture  point,  noting  that  the  surface 
curvature,  k,  takes  a value  of  one  on  the  spherical  part 
and  a value  of  zero  on  the  conical  part  and  also  that 
the  flow  properties  are  continuous  through  this  juncture 
point.  This  procedure  is  similar  to  that  adopted  for 
the  inviscid  set  of  equations.  Note  also  that  the  jump 
conditions  associated  with  the  shock  shape  derivatives 
would  remain  the  same  as  those  for  the  inviscid  case. 

Thus  a proper  physical  behavior  of  the  full  viscous  shock 
layer  equations  at  the  sphere/cone  juncture  point  is 
summarized  as  follows: 

1.  The  flow  variables  are  continuous  at  the  sphere/ 
cone  juncture  point. 

2.  The  use  of  a surface  coordinate  system  intro- 
duces discontinuities  in  the  flow  gradients 
relative  to  surface  distance  everywhere  across 
the  shock  layer  at  the  juncture  point. 

3.  Independent  of  the  choice  of  the  coordinate 
system,  inviscid  theory  predicts  discontinuities 
in  flow  gradients  only  at  the  surface  at  the 
sphere/cone  juncture  point.  However  the  viscous 
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flow  analysis  of  Messiter  [7]  indicates  that 
in  the  limiting  case  of  very  high  Reynolds 
number,  this  discontinuity  would  be  smoothed 
out  by  the  sublayer  interaction  effect  within 
the  inner  scale  length. 

4.  .Within  the  viscous  layer  the  gradient  discontinui- 
ties due  to  the  choice  of  the  coordinate  system 
would  tend  to  drop  out  of  the  lead  order  viscous 
equations  as  the  Reynolds  number  tends  to  infinity. 

3 . Thin  Layer  Analysis 

Many  studies  [13,  14]  in  the  past  have  used  the  thin 
layer  version  of  the  full  shock  layer  equations  to  predict 
flow  properties  within  the  shock  layer  region  for  analytic 
bodies  such  as  spheres,  paraboloids  and  hyperboloids  at 
high  Mach  number.  The  simplifying  assumptions  inherent 
in  the  thin  shock  layer  approximations  cause  a change  in 
the  character  of  the  governing  equations  and  thus  of  the 
juncture  point  analysis  presented  above.  In  order  to 
analyze  this  set  of  equations,  the  inviscid  form  of  these 
equations  are  first  considered  here.  Attention  is  first 
drawn  to  the  characteristics  of  these  equations.  In  the 
surface  coordinate  system  the  inviscid  thin  shock  layer 
equations  are  given  by  equations  (17,  18,  20,  21)  and  the 
normal  momentum  equation  is  given  as 
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Pn(l+Kn)  - KpU2  = 0 (31) 

These  sets  of  equations  can  be  shown  to  be  parabolic  in 
nature  indicating  that  the  characteristics  of  the  flow  are 
perpendicular  to  the  surface  of  the  body  (see  Appendix  C) . 
This  suggests  that  for  the  thin  shock  layer  equations 
information  from  the  body  surface  is  propagated  along  a 
line  perpendicular  to  the  body  surface  unlike  the  full 
shock  layer  equations  where  the  characteristic  lines  are 
inclined  at  the  local  Mach  angle  of  the  flow.  For  this 
reason  it  is  obvious  that  any  discontinuity  in  the  flow 
derivatives  or  otherwise  discontinuity  at  the  sphere/cone 
juncture  would  be  felt  all  across  the  shock  layer 
immediately  at  the  juncture  point.  This  physical  behavior 
of  the  inviscid  thin  shock  layer  equations  is  significantly 
different  from  their  full  shock  layer  counterpart  and  would 
be  expected  to  manifest  itself  rather  dramatically  in 
the  solutions  obtained. 

In  order  to  further  study  the  behavior  of  the  thin 
shock  layer  equations  at  the  sphere/cone  juncture  point 


* It  can  also  be  shown  from  the  analysis  of  the  characteris- 
tics of  the  full  shock  layer  equations  that  their  charac- 
teristics tend  to  become  perpendicular  to  the  body  surface 
in  the  limit  as  y -*■  1,  i.e.  the  thin  shock  layer 
approximation  is  approached  (see  Appendix  C) . 
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attention  is  now  directed  to  the  normal  momentum  equa- 
tion given  above.  This  shows  that  the  pressure  will 
be  a constant  across  the  shock  layer  on  the  conical 
portion  of  the  body  whereas  a pressure  variation  will 
be  encountered  on  the  spherical  portion.  This  can 
only  occur  if  a jump  is  allowed  in  the  pressure  at  the 
sphere/cone  tangency  point.  It  is  to  be  noted  here 
that  this  need  for  a discontinuity  in  pressure  is  valid 
for  inviscid  as  well  as  viscous  flows  since  the  normal 
momentum  equation  under  the  thin  layer  approximations 
remains  unaltered.  However  it  is  of  interest  to  note 
that  since  under  the  inviscid  thin  layer  approximation 
the  information  is  propagated  normal  to  the  body 
surface,  the  bow  shock  wave  is  expected  to  feel  the 
presence  of  the  sphere/cone  juncture  point  and  its 
attendant  pressure  discontinuity  immediately  above  the 
juncture  location. 

It  is,  thus,  important  to  the  whole  structure  of 
the  flow  that  the  nature  of  the  thin  layer  solutions  be 
delineated.  To  determine  the  thin  layer  "jump  condition", 
the  approach  taken  here  is  to  revert  to  the  integral  form 
of  the  full  governing  equations,  and  to  assess  the 
influence  of  the  thin  layer  approximations  on  the 
generalized  jump  conditions  so  obtained.  To  do  this  one 
must  first  identify  the  "thin"  layer  terms  in  the  general 
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analysis,  identify  their  contributions  to  the  integral 
formation  of  the  governing  equation,  and  then  make  the  thin 
layer  assumption.  In  a manner  similar  to  the  full  shock 
layer  equations,  an  element  of  infinitisimally  small  size 
is  considered  in  this  coordinate  system  as  shown  in 
Appendix  (A) . The  integral  form  of  the  momentum  equation 
(All)  when  evaluated  for  this  element,  is  shown  to  be 
equation  (A12)  when  no  approximation  is  made.  It  is  now 
necessary  to  neglect  from  this  equation  those  terms  which 
lead  to  thin  shock  layer  approximations.  However,  at 
this  point  those  terms  which  are  neglected  in  making 
these  assumptions  are  unknown.  Therefore,  one  must 
extend  this  derivation  to  obtain  the  differential  form 
of  the  governing  equations  and  track  back  those  terms 
which  are  neglected  when  thin  shock  layer  assumptions  are 
made.  This  is  achieved  in  Appendix  (A) . The  equations, 
so  obtained,  are  given  as: 

Pl  + Pl  u*  = p2  + p2  (32) 


P1  U1  “ p 2 U2 

2 2 

U1  u2 

hl  + 2 = h2  + 2 


(33) 

(34) 
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Note  that  these  equations  differ  from  the  corresponding 
equations  for  the  full  shock  layer  form  only  in  the  one 
respect  that  the  v component  of  velocity  does  not  appear 
in  the  present  form  and  thus  is  unrestricted  in  its 
jump  behavior  across  the  juncture  point. 

The  equations  (32-34)  present  the  set  of  conditions 
that  must  be  met  at  the  sphere/cone  juncture  by  all  the 
flow  variables  in  order  to  accommodate  the  pressure  jump 
predicted  by  the  normal  momentum  equation.  These  equations 
are  quite  similar  to  the  "shock  discontinuities"  of  a 
perfect  gas  except  that  they  do  not  contain  the  v-component 
of  velocity.  The  fact  the  v-component  of  velocity  does 
not  appear  in  these  jump  conditions  is  not  surprising 
since  a thin  shock  layer  approximation  tacitly  assumes  the 
v-component  of  velocity  to  be  small  compared  to  the 
longitudinal  velocity  and  as  such  plays  a secondary  role 
in  the  conservation  laws. 

The  admissibility  of  jump  conditions  in  the  flow 
variables  at  the  juncture  point  renders  the  physics  of 
the  flow  rather  complex  in  this  region.  At  the  bow  shock, 
such  jump  conditions  would  be  expected  to  produce  a 
discontinuous  shock  slope  at  the  juncture  point.  It 
would  appear  that  the  juncture  discontinuity  propagates 
normal  to  the  body  surface  through  the  local  characteristics 
directly  to  the  bow  shock  shape. 
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It  is  not  really  clear  how  the  discontinuities  in 
the  flow  variables  themselves  are  accommodated  by  the 
thin  shock  layer  equations.  Surely  discontinuities  in  the 
gradients  are  to  be  expected  but  more  confusing  is  the 
anticipated  behavior  of  the  viscous  flow  regions.  Here 
it  is  not  clear  whether  discontinuous  solutions  can  exist 
since  preliminary  analysis  would  seem  to  indicate  that 
they  would  be  of  the  subsonic  "expansion  shock"  type. 
Further  study  of  this  point  is  warranted  but  shall  be 
deferred  from  the  current  effort. 
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IV.  NUMERICAL  METHOD 

1.  General  Considerations 

Several  methods  have  been  presented  for  solving  the 
"thin"  shock  layer  version  of  the  more  general  viscous 
shock  layer  equations  [13,  14].  These  approaches  have 
two  limitations.  First  they  are  based  on  the  assumption 
that  the  pressure  gradient  normal  to  the  body  surface  is 
established  entirely  by  centrifugal  effects,  and  second, 
that  the  shock  wave  lies  parallel  to  the  body  surface. 

In  an  attempt  to  remove  these  limitations,  methods  have 
been  developed  [2,  15]  for  addressing  the  full  shock 
layer  equations  through  a relaxation  process  wherein  the 
thin  shock  layer  assumptions  are  removed  by  an  iterative 
process.  While,  in  general,  such  methods  have  been 
successful,  they  encounter  difficulty  whenever  the  shock 
layer  thickness  becomes  large.  This  difficulty  usually 
manifests  itself  as  a divergent  behavior  in  the  iteration 
scheme.  In  an  attempt  to  overcome  this  problem,  a new 
relaxation  scheme  was  developed  in  Reference  [16]  where 
an  initial  solution  was  relaxed  in  an  artificial  time  like 
manner  toward  the  sought  after  "steady  state"  solution. 

In  this  sense  the  approach  is  similar  to  the  relaxation 
scheme  presented  by  Davis  [2],  Davis  and  Nei  [17]  and 
Srivastava,  Werle  and  Davis  [18],  the  primary  difference 
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being  the  manner  in  which  the  "new  guess"  on  the  solution 
is  defined  after  a given  step  in  the  relaxation  process. 
The  method  so  developed  has  been  found  to  work  well  for 
analytic  bodies  such  as  paraboloids,  hyperboloids  and 
spheres . 

Application  of  this  approach  to  nonanalytic  bodies 
such  as  spherically  blunted  cones  encounters  numerical 
difficulties  at  a sphere/cone  juncture  point  where  the 
longitudinal  flow  derivatives  in  a surface  coordinate 
system  undergo  discontinuous  changes.  The  method  of 
solution  presented  here  represents  an  adaptation  of  the 
earlier  time  like  relaxation  scheme  [16]  to  problems 
with  imbedded  discontinuities  in  the  flow  derivatives. 

In  order  to  demonstrate  the  present  approach,  the 
s-momentum  equation  of  the  viscous  shock  layer  equations 
is  first  rewritten  in  the  form 


32s 


3n' 


+ So  + 3.,  P*  + g 

*1  3n  2 dg2  3 ds  *4 


s = 0 
P5  35  ° 


(35) 


where  3^#  32?  B3/  84  and  65  can  be  obtained  from  Reference 
[16]  and  are  given  in  Appendix  (D) . 

The  present  time  like  relaxation  scheme  utilizes  a 
two  step  process  in  which  the  first  step  of  the  method  is 
somewhat  similar  to  an  alternating  direction  implicit 
method  and  it  yields  the  flow  variables  in  the  shock  layer 


39 


AEDC -TR-77-20 


region  while  the  second  step  is  used  to  update  the  shock 
shape  itself.  This  scheme  can  be  demonstrated  through 
the  s-momentum  equation  (7)  written  in  a two  step  time 
formulation  (see  Fig-.  3)  as. 


First  Sweep  * 


.2— * * — * . 2_n 

3 u-  + R*  *,3  R 

2 + ei  3n  B2l._2 


an 


Second  Sweep  n+1 


8s 


if -o 

(36a) 


* a2Rn+1 
'2 


* *«n+l 


_ 8 + 8*  8K 

S2  8 1 - + S3  8s 


n+1 


+ «5  If  + »4>*  = 0 


(36b) 


Note  that  the  "steady  state"  version  of  these  equations  are 

i 

precisely  the  "full"  shock  layer  equations. 

The  boundary  conditions  associated  with  star  sweep 
equation  (36a)  are  typical  no  slip  conditions  (2a,  2b)  at 
the  surface  and  Rankine-Hugoniot  conditions  (3a-3g)  at  the 
shock  location.  However,  the  boundary  conditions  associated 
with  the  final  sweep  are  the  same  as  those  used  in 
Reference  [16]  and  are  given  as. 


i)  At  s = 0 R = 0 


(37a) 


ii)  At  s 


smax 


Rn+1  = R* 

max  max 


(37b) 
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There  are  two  points  of  interest  produced  by  the 
sphere/cone  curvature  discontinuity.  First  note  that  the 
coefficients  in  the  s-momentum  equation  (7)  contain  the 
surface  curvature,  which  itself  undergoes  a discontinuous 
change  from  a value  of  one  on  the  spherical  body  to  a 
value  of  zero  on  the  conical  body  at  the  sphere/cone 
juncture  point.  This  then  causes  the  B coefficients  of 
equation  (36)  to  undergo  discontinuous  changes  at  the 
juncture  point.  As  a second  point  of  interest,  it  is 
also  noted  that  the  derivatives  of  the  shock  shape  with 
respect  to  s,  explicitly  appear  in  the  governing  equations, 
such  as  in  equation  (36) . These  shock  derivatives  have 
been  shown  to  undergo  discontinuities  at  the  juncture 
point  in  Section  III.  This  jump  condition  on  the  first 
derivative  of  the  shock  shape  is  given  by  equation  (26) . 

A similar  jump  condition  on  the  second  derivative  of  the 
shock  shape  can  be  obtained  through  use  of  the  governing 
differential  equations. 

As  a result  then  it  is  seen  that  the  governing 
equation  (36)  contains  both  flow  coefficients  and  shock 
derivatives  which  undergo  discontinuous  change  at  the 
juncture  point  which  necessarily  produce  solutions  with 
discontinuous  gradients  (see  Section  III) . Such  results 
obviously  require  modifications  in  order  to  obtain 
numerical  solutions  of  this  set  of  governing  equations 
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to  assure  that  finite  differencing  is  not  done  across 
such  discontinuous  regions.  Note  that  the  numerical 
difficulties  associated  with  equation  (36a)  can  be  over-  ■ 
come  by  structuring  the  finite  difference  grid  system 
with  a point  at  the  juncture  thereby  avoiding  differencing 
across  the  discontinuity.  However  difficulty  is  still 
encountered  in  the  second  step  of  the  solution  process 
due  to  the  discontinuities  that  occur  in  the  shock  shape 
derivatives.  The  occurrence  of  these  discontinuities 
requires  the  use  of  special  difference  relations  at  the 
juncture  point.  It  will  be  shown  here,  through  a model 
problem  representing  the  second  step  of  the  present 
numerical  scheme  that  this  difficulty  can  be  overcome  if 
the  difference  form  of  the  differentials  are  formulated 
such  that  they  comprehend  the  juncture  jump  conditions. 

2 . A Model  Problem 

The  governing  equation  for  a model  problem  and  its 
associated  boundary  and  jump  conditions  are  formulated  here 
in  order  to  demonstrate  the  concepts  associated  with  the 
viscous  shock  layer  solution  for  spherically  blunted  cones. 

The  governing  equation  for  this  model  problem  is  taken 
to  be  analogous  to  the  second  step  of  the  viscous  shock 
layer  scheme  and  is  given  as 
2 

+ “i  it  + “2  R + “3  = 0 <38> 
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Note  that  the  coefficients  a^,  a 2'  “3  in  the  above 
equation  are  to  be  selected  so  that  this  equation  models 
the  second  step  of  the  viscous  shock  layer  scheme.  This 
would  require  that  either  some  or  all  of  these  coeffi- 
cients undergo  a discontinuous  change  at  the  model  juncture 
point  in  the  solution  region.  For  present  purposes  the 
coefficients  are  taken  to  be  one  set  of  constants  in  the 
region  ahead  of  the  juncture  and  a different  set  of 
constants  aft  of  the  juncture.  A comparison  with  the 
second  step  of  the  viscous  shock  layer  scheme  shows  that 
the'  coefficient  a2  does  not  encounter  any  jump  whereas 
and  a3  do  encounter  jumps  in  their  magnitudes  at  the 
sphere/cone  tangency  point.  Thus,  the  present  model 
problem  is  set  up  such  that  the  coefficients  c^,  a2/  a3 
take  constant  values  (corresponding  to  the  spherical  section) 
in  one  region  and  different  constant  values  (corresponding 
to  the  conical  section)  aft  of  the  juncture  location. 

The  boundary  conditions  to  be  applied  to  this  model 
problem  are  established  to  closely  correspond  to  the 
second  step  of  the  viscous  shock  layer  solution.  These 
then  are  given  as 

R = 0 at  s = 0 (39a) 

R = Ri  at  s " smax  <39b> 
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In  addition  to  these,  jump  conditions  analogous  to  those 
discussed  for  the  viscous  shock  layer  equations  must  be 
established  and  applied  at  the  juncture  point. 

At  the  juncture  it  is  required  that  the  function 
R,  be  continuous  but  that  the  first  derivative,  dR/ds, 
be  represented  by  the  relation 


<s> 


-s. 


= Ki  <af> 


(40) 


+s_ 


Dump  "jump 

Note  that  the  value  of  K1  will  be  determined  here  to  match 

the  jump  actually  encountered  by  the  viscous  shock  layer 

shock  derivative  at  the  sphere/cone  tangency  point. 

From  the  governing  equation  (38) , it  is  found  that  the 
. . 2 2 

second  derivative,  d R/ds  , also  undergoes  a jump  at  the 
juncture  point  as  given  by  the  relation 


<^S) 


ds  -s 


- (a!S) 


Dump 


,,  + K,<5§) 

ds2  +s.  _ 2 as  +8 

Dump 


+ K 


jump 


3 

(41a) 


where 


K2  = (al>+s  “ Kl(ai>-S 

j ump  j ump 


K3  = ^^-s.  " ^a3^  +s  . 

Dump  ®]ump 


(41b) 


The  exact  solution  for  the  above  equations,  associated 
boundary  conditions  and  jump  requirements  can  be  found  easily 
and  is  given  as 


R = A ems  + B e“ns  - a3/a2 
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for  either  region,  where 

“5 

m = - a^/2  + /aj/4  - c*2  >0  for  c*2  < 0 

n = cij/2  + /o£/4  - a2  >0  for  a2  < 0 

* 

This  gives, 

R = (<*3/°^)  ems  + B[e  nS  - ems]  - a3^a2  (*2) 


when 

Also, 


0 < s < -s 


jump 


D, (s-s  ) D,  s-2s^  ) -D, s 

„ « 1 max'  _ r 1 max  _ 1 

R = R^  e - B^le  - e 

D, (s-s  ) 

+ a2/a2  [e  - 1J 


(43) 


for 


+s 


jump  — 


s 


s_ 


max 


where  for  = 0 in  this  region.  Equations 

(42,  43)  contain  two  undetermined  constants  B and  which 
can  be  determined  by  using  the  jump  condition  (40)  and  the 
condition  that  R be  continuous  at  the  juncture  point. 


3.  Model  Problem  Finite  Difference  Formulation 

The  finite  difference  form  of  the  differentials  in 
equation  (38)  requires  special  attention  at  the  juncture 
point  in  order  to  account  for  the  jump  conditions  associated 
with  the  various  derivatives.  This  formulation  is  shown 
through  the  figure  below  where  typical  mesh  points  are 
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shown  with  a jump  occurring  immediately  ahead  of  point  3. 

Juncture  Point 


12  3 4 


Mesh  point  "a"  is  located  immediately  ahead  of  the  juncture 
point.  While  formulating  the  difference  form  of  the  deri- 
vatives at  point  "3" , a Taylor's  series  representation 
which  avoids  any  series  expansion  across  the  discontinuity 
is  utilized.  For  this  case,  then,  a Taylor'B  series  expan- 
sion is  used  from  point  3 to  4 and  from  point  "a"  to  2. 
Points  "a"  and  3 across  the  discontinuity  are  related 
through  the  jump  conditions  (40,  41).  This  procedure  yields 
after  proper  manipulation  (see  Appendix  E)  the  following 
forms  of  the  difference  representations  of  the'  differentials 


in  equation  (38)  at  the  juncture  point. 

AK, 


,dR*  _ 

W " 


R4  ~ R2 


3 A (1+K1~  f K2)  2(1+Kx-  | K2) 


+ 0(A2) 


(44a) 


2K- 


(d2R)  o (R2~R3)  + (K1~  2 K2*  (R4"R3)  _ “'3 

dl7  3 (Ai/2X1+K1-  | K2)  A2  (1+Ki-|  K2) 


+ 0 (A  ) 


(44b) 


Note  is  made  here  that  these  expressions  contain  the  constants 
associated  with  the  jump  conditions  (40,  41)  and  that  they 
reduce  to  their  proper  central  difference  representation  for 
zero  strength  jump  (e.g.,  K^l,  K2=0  and  K3=0)  . A formally 
second  order  accurate  representation  can  similarly  be 
formulated  for  the  second  derivative  by  evaluating  the 
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error  term  of  equation  (44b)  using  the  governing  differen- 
tial equation  at  the  point  where  the  jump  is  taking  place 
in  the  solution  regime.  The  error  term  in  Equation  (44b) 
is  given  as 


E = 


4 . (K1-  I 

3 1 


k2) 


- R 


d+Ki  - 7 


k2) 


-]  + 


0 (A2) 


For  the  model  problem  the  first  term  (order  A)  of  the  error 
expression  can  be  evaluated  by  first  differentiating  the 
governing  differential  equation  and  then  evaluating  it  on 
the  two  sides  of  the  point  where  the  jump  in  derivatives 
occurs.  This  yeilds 


Ra"  + “la  Ra'  + “2a  = 0 


(45a) 


and 


f I V 


R3  + “l3  R3  + “23  R3  ‘ 0 


(45b) 


Manipulation  of  equations  (45a, b)  results  in 


I 9 I 


(K1  “ 7 K2)R3  “ - £ Kn)a,  R 


1 7 *2,al. 


+ «1.  E»’  + “2.  K - <K1  - I K2>  “2.^2 


(46! 


A substitution  of  this  expression  in  (44b)  would  result  in 
a second  order  accurate  second  derivative  at  the  juncture 
point.  For  further  details  see  Appendix  (E) . 
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For  the  sake  of  identification  in  the  rest  of  this 
text,  the  first  of  these  above  difference  expressions  (44) 
will  be  referred  to  as  the  "first  order  accurate  scheme", 
even  though  only  the  second  derivative  at  the  juncture 
point  experiences  a formally  first  order  error.  The 
second  of  these  where  the  second  derivative  is  also 
formally  second  order  accurate  will  be  referred  to, 
here,  as  the  "second  order  accurate"  scheme. 

Before  the  numerical  scheme  outlined  above  is  used 
in  the  viscous  shock  layer  equations,  a test  of  this 
scheme's  ability  to  approximate  the  exact  solution  is 
essential. 

4 . Model  Problem  Numerical  Results 

Figure  4 shows  the  variable  "R"  as  a function  of 
distance  for  the  model  problem.  The  coefficients  c^, 
a2'  a3  the  mode^  equation  were  chosen  to  approxi- 
mately represent  the  sphere/cone  juncture  point  of  a 
40°  half  angle  spherically  blunted  cone.  Severe  error 
is  seen  in  the  numerical  solution  when  the  jump  effects 
are  ignored.  However,  when  proper  jump  effects  are 
accounted  for,  the  exact  solution  is  virtually  recovered 
by  the  present  numerical  scheme.  Figure  5 shows  the  first 
derivative  for  the  same  problem.  It  is  noted  that  the 
discontinuity  predicted  by  the  exact  solution  is  virtually 
captured  exactly  by  the  numerical  scheme  if  proper  jump 
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effects  are  accounted  for  in  the  numerical  scheme.  Figure 
6 shows  the  second  derivative,  d R/ds  , as  a function  of 
surface  distance.  It  is  seen  also  from  this  figure  that 
large  numerical  errors  are  present  when  the  effects  of 
the  jump  are  ignored.  However,  the  exact  solution  is 
accurately  recovered  when  proper  jump  effects  are 
included. 

For  the  sake  of  completeness,  the  model  problem  was 

also  solved  using  the  "second  order  accurate  scheme"  of 

equation  (44)  . It  was  found  that  for  this  case  the 

numerical  difference  could  not  be  detected  to  the  scale 

of  plot  shown  in  Figures  4 through  6.  In  order  to  clarify 

this  point,  further  studies  were  undertaken.  Figure  7 

2 2 

shows  the  model  problem  shock  curvature,  d R/ds  , at 
the  junction  point  as  a function  of  the  step  size  "As" 
for  the  "first  order"  and  "second  order"  numerical  schemes. 
Note  here  is  made  of  the  fact  that  the  so-called  "first 
order  scheme"  employs  a formally  first  order  accurate 
representation  of  the  second  derivative  at  only  one  point 
in  the  mesh  system,  i.e.  at  the  juncture  point.  It  is  seen 
from  Figure  7 that  the  "second  order  scheme"  shows  a 
parabolic  behavior  as  the  step  size  "As"  is  reduced,  as 
one  would  expect.  However,  the  "first  order  scheme"  does 
not  show  a linear  dependence  on  As  but  rather  a parabolic 
behavior.  To  further  detail  these  results,  Figure  8 replots 
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these  curves  against  the  square  of  the  step  size  "As". 

This  figure  clearly  indicates  that  both  the  "first"  and 
"second"  order  schemes  yield  results  that  approach  the 
exact  solution  as  though  they  were  second  order  accurate. 
The  explanation  of  these  results  is  given  in  detail  in 
Appendix  (F)  where  it  is  verified  that  it  is  basically  a 
manifestation  of  the  fact  that  a local  truncation  error 
of  order  As  at  a finite  number  of  points  in  a finite 
difference  mesh  does  not  necessarily  produce  a first  order 
global  error. 

It  is,  therefore,  established  that  the  so-called 
"first  order  accurate  scheme"  is  essentially  second  order 
accurate.  For  this  reason  it  was  found  unnecessary  to  use 
the  "second  order  accurate  scheme"  for  present  purposes. 

5.  Application  to  the  Full  Viscous  Shock  Layer  Equations 

The  finite  difference  expressions  so  developed  can 
now  be  applied  to  the  solution  of  the  full  viscous  shock 
layer  equations  for  hypersonic  flow  past  spherically 
blunted  cones.  A detailed  description  of  the  method  for 
evaluating  the  jump  conditions  associated  with  the  shock 
wave  derivatives  at  the  sphere/cone  juncture  point  for  the 
full  viscous  shock  layer  equations  is  presented  in 
Appendix  (G) . It  is  shown  here,  that  the  proper  jump 
condition  associated  with  the  first  derivative,  dR/ds, 
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at  the  juncture  point  can  be  estimated  by  purely  geometrical 
considerations  and  the  fact  that  the  shock  wave  shape,  R, 
itself  is  smooth  at  this  point.  However,  the  jump  condition 
associated  with  the  second  derivative,  d R/ds  , must  be 
obtained  by  using  the  differential  equation  itself  on  the 
two  sides  of  the  sphere/cone  tangency  point  in  a manner 
similar  to  that  adopted  for  the  model  problem  discussed 
earlier.  This  is  shown  further  in  Appendix  (G) . It  is 
shown  in  this  appendix  that  care  must  be  used  while 
evaluating  the  jump  conditions  in  the  viscous  shock  layer 
code.  For  present  purposes,  the  jump  conditions  were 
evaluated  at  the  first  mesh  point  away  from  the  wall. 

6 . Overall  Method  of  Solution 

The  overall  method  of  solution  for  the  full  viscous 
shock  layer  equations  is  as  follows.  An  initial  guess 
was  first  made  on  the  shock  shape.  Based  on  this  guess 
the  first  and  second  derivatives  of  the  shock  distance 
were  computed  using  central  differences  at  points  away 
from  the  sphere/cone  tangency  point.  However,  since 
jump  conditions  on  the  shock  derivatives  are  not  known 
initially  at  the  sphere/cone  juncture  point,  a second  order 
accurate  three  point  backward  difference  schemes  was  used 
on  the  spherical  part  and  a three  point  forward  difference 
scheme  was  used  on  the  conical  part  for  the  first 
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derivative  of  the  shock  stand  off  distance  in  order  to 
avoid  any  differencing  across  the  juncture  point.  In  a 
like  manner,  four  point  second  order  accurate  schemes 
were  used  for  the  second  derivative  of  the  shock  stand 
off  distance  at  the  juncture  point.  The  star  sweep  equations 
were  then  solved  by  starting  at  the  stagnation  point, 
where  both  3u/3f;  and  3t/3f;  vanish, thus  reducing  the 
governing  equations  to  ordinary  differential  equations. 

The  first  equation  solved  was  the  energy  equation  so 
that  thereafter  all  quantities  such  as  viscosity  related 
to  temperature  could  be  evaluated.  Next,  the  s-momentum  . 
equation  was  integrated  to  determine  a u-velocity  profile, 
and  then  the  continuity  equation  was  solved  to  determine 
first  the  shock  stand-off  distance  from  equation  (10) 
and  then  the  v-component  of  the  velocity  from  equation 
(16) . Finally  equation  (11)  was  integrated  to  determine 
the  local  pressure  level.  The  coefficients  in  the 
governing  equations  were  then  reevaluated  using  the  new 
flow  variables.  Repetition  of  the  above  steps  at  a 
given  station  continued  until  the  solution  converged. 

The  method  then  stepped  along  the  body  surface  and 
iterated  at  each  station  to  achieve  converged  solutions. 

To  accelerate  the  convergence  process,  the  previous 
station  values  of  the  profiles  were  used  at  each  new  step 
as  a first  guess.  One  difficulty  encountered  during  this 
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iteration  scheme  was  the  presence  of  an  oscillatory 
behavior  of  the  normal  velocity  component,  v,  at  some 
station  in  the  s-direction  [Ref.  16] . This  oscillatory 
behavior  of  the  physical  quantities  was  overcome  by  an 
under-relaxation  scheme  as  shown: 

w = F^  w^  + (1-F1)w2 

where  w^  is  the  most  recently  calculated  physical  quantity 

and  w2  is  the  value  obtained  from  the  previous  calculated 

value  of  this  quantity.  It  was  observed  that  a value 

of  F^  of  0.2  to  0.4  produced  convergence  in  most  cases 

considered.  In  general  it  was  also  found  that  such  an 

under-relaxation  technique  was  needed  only  for  the  pressure 

and  v-component  of  velocity. 

Once  the  above  method  had  passed  over  the  entire  mesh 

the,  second  sweep  equations  were  invoked.  The  final  sweep 

equation  (36b)  was  then  solved  using  the  two  boundary 

conditions  of  equations  (37a,b) . No  iteration  of  the 

final  sweep  equation  was  required  since  it  is  linear. 

However  note  that  the  final  sweep  equation  requires  the 

necessary  jump  conditions  associated  with  the  first 

derivative  of  the  shock  stand  off  distance,  dR/ds  and 

2 2 

also  that  associated  with  the  second  derivative,  d R/ds  . 
These  jump  conditions  were  evaluated  using  the  flow 
properties  obtained  in  the  star  sweep  calculations.  The 
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final  sweep  equation  was  required  since  it  is  linear. 
However  note  that  the  final  sweep  equation  requires  the 
necessary  jump  conditions  associated  with  the  first 
derivative  of  the  shock  stand  off  distance,  dR/ds  and 
also  that  associated  with  the  second  derivative,  d R/ds  . 
These  jump  conditions  were  evaluated  using  the  flow 
properties  obtained  in  the  star  sweep  calculations.  The 
shock  shape  obtained  from  the  final  sweep  was  used  then 
to  solve  the  next  star  sweep  in  time.  The  procedure 
continued  in  time  until  two  alternate  final  sweeps 
converged  to  a desired  degree  of  accuracy.  Appendix  (H) 
discusses  further  details  of  the  computer  program  used 
to  obtain  the  present  numerical  results. 

7.  Grid  Sizes  for  Shock  Layer  Solution 

The  following  normal  step  sizes  distributions  were 

t 

used  in  the  finite  difference  solution  of  the  full  viscous 
shock  layer  equations  for  the  cases  presented  in  the 
following  section. 


Re  = 1.515 
00 

x 103 

n range 

An 

0.0  - 0.050 

0.001 

0.05  - 0.65 

0.015 
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00 

105 
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0.0  - 0.005 

0.0001 

0.005  - 0.50 

0.005 

0.50  - 1.0 

0.0099 

54 


AEDC-TR-77-20 


V.  RESULTS  AND  DISCUSSION 

The  general  analysis  and  the  numerical  techniques 
discussed  earlier  were  used  to  obtain  the  solutions  of 
the  full  viscous  shock  layer  equations  for  hypersonic 
flow  past  various  spherically  blunted  cones  in  order  to 
test  the  reliability  of  this  technique.  Since  the  interest 
in  the  present  study  was  centered  on  the  sphere/ cone 
tangency  region,  numerical  solutions  were  generated  only 
to  about  2-3  nose  radii  downstream  of  the  stagnation 
point  for  a range  of  large  as  well  as  small  cone  angles 
of  the  spherically  blunted  cones.  Numerical  solutions 
were  obtained  for  a wide  range  of  cone  half  angles  from 
30°  to  0°  at  various  test  conditions  corresponding  to 
available  data  and  other  calculations. 

Figure  9 shows  the  surface  pressure  distribution 
obtained  here  for  a 30°  half  cone  angle  spherically  blunted 
cone  at  a free  stream  Mach  number,  = 10,  free  stream 
Reynolds  number,  Re^  = 3 x 10 ^ and  a wall  to  stagnation 
temperature  ratio,  = 0.05.  These  test  conditions 

were  chosen  in  order  to  compare  the  predicted  numerical 
results  with  the  inviscid  solutions  of  Inouye  et  al. 

[8]  for  the  same  body.  The  present  calculations  for  this 
case  were  made  using  a variable  normal  step  size,  An, 
which  ensured  at  least  10-15  mesh  points  within  the 
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boundary  layer  regime  while  the  longitudinal  step  size, 
is,  was  selected  such  that  a mesh  point  of  the  numerical 
scheme  coincided  with  the  sphere/cone  juncture  point. 

The  time  step  size,  it,  was  taken  to  be  3.5. 

Figure  9 shows  that  when  the  shock  jump  conditions 

are  accounted  for  in  the  finite  difference  formulation 

as  shown  earlier,  the  predicted  surface  pressure  compares 

well  with  the  inviscid  solution.  It  is  of  interest  to 

note  here  that  the  discontinuity  in  the  surface  pressure 

gradient  at  the  sphere/cone  juncture  point  predicted  by 

the  inviscid  theory  is  virtually  reproduced  by  the 

present  viscous  model  for  this  very  high  Reynolds  number 

case.  However,  in  the  light  of  the  analysis  of  Messiter 

and  Hu  [7]  for  a simple  two-dimensional  flow  with  a 

curvature  discontinuity,  one  might  have  anticipated  a 

viscous  smoothing  of  the  discontinuity  in  the  surface 

pressure  gradient  at  the  juncture  point.  It  is  now  clear 

* 

from  the  discussion  of  Section  III  that  the  viscous 
smoothing  for  this  problem  is  of  a very  mild  nature  and 
occurs  over  such  a short  distance  that  it  is  not  seen 
to  the  scale  of  the  present  calculations. 

Figure  9 also  presents  two  other  numerical  results 
for  the  same  test  conditions.  One  of  these  is  the  case 
where  a numerical  solution  of  the  full  viscous  shock 
layer  equations  was  obtained  by  ignoring  the  relevant 
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jump  conditions  associated  with  the  shock  derivatives  at 
the  juncture  point  and  by  adjusting  the  mesh  of  the  nu- 
merical scheme  such  that  the  juncture  point  lies  between 
two  mesh  points.  This  case  then  allows  an  assessment 
of  the  numerical  errors  that  are  introduced  in  the 
computational  method  when  one  ignores  the  relevant  jump 
effects  associated  with  the  juncture  point.  Note  that 
the  numerical  errors  are  large  in  the  juncture  region 
and  tend  to  diminish  away  from  the  juncture  region.  The 
second  case  shown  in  Figure  9 is  similar  to  the  first 
except  that  the  mesh  system  was  aligned  so  that  a mesh 
point  coincided  with  the  sphere/cone  tangency  point. 

This  case  allows  an  assessment  of  the  importance  of  the 
shock  jump  conditions  on  the  surface  properties. 

Apparently  these  are  of  a dominant  nature  in  this  region 
of  the  flow. 

Figure  10  shows  the  surface  heat  transfer  distribution 
for  this  case  under  identical  flow  conditions.  It  is 
observed  that  the  erratic  behavior  of  the  computational 
results  persist  when  the  jump  effects  are  ignored  whereas 
the  inclusion  of  these  effects  tend  to  eliminate  this 
erratic  behavior  completely.  From  equation  (26)  it  is 
seen  that  the  jump  conditions  associated  with  the  shock 
derivatives  tend  to  increase  in  magnitude  as  the  cone 
angle  for  the  spherically  blunted  cones  is  reduced.  It 
is,  therefore,  pertinent  to  test  this  computational 
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technique  for  lower  cone  angles  in  order  to  establish 
the  generality  of  this  scheme.  Further  numerical  solu- 
tions were  obtained  for  lower  cone  angles  ranging  from 
20°  to  0°  in  order  to  assess  the  influence  of  the  jump 
effects  on  the  surface  properties.  Difficulties  were 
encountered  here  while  attempting  to  reduce  the  cone 
angle  mainly  due  to  the  choice  of  the  initial  shock  shape 
for  such  bodies.  This  was  overcome  here  by  reducing 
the  cone  angle  in  increments  of  about  5°  with  the  number 
of  mesh  points  between  the  juncture  and  stagnation  point 
kept  fixed.  This  resulted  in  an  increase  in  longitudinal 
step  size,  As,  as  the  cone  angle  was  reduced  but  this 
technique  was  found  to  work  well  for  all  cases  that  are 
presented  here.  It  should  also  be  noted  that  care 
was  exercised  to  at  least  include  10-15  points  within  the 
boundary  layer  while  selecting  the  normal  step  size,  An. 
Figure  11  through  15  present  the  results  for  such  a 
calculation  for  cone  angles  ranging  from  20°  to  0°.  It 
is  seen  from  Figures  11  and  13,  which  show  the  surface 
pressure  distributions, that  the  results  of  such  a cal- 
culation compare  well  with  the  inviscid  solution,  when 
the  proper  jump  conditions  associated  with  the  shock 
derivatives  are  included  in  the  solution  scheme.  These 
figures  also  show  the  case  when  such  jump  effects  are 
ignored  in  the  calculation,  thereby  causing  rather 
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erratic  behavior.  Figures  12  and  14  show  similar  results 
for  the  surface  heat  transfer  rates.  Note  that  Figure  15 
for  a 0°  cone  angle  (spherical-cylinder  body)  does  not 
contain  the  numerical  results  corresponding  to  the  no 
jump  case.  This  is  because  the  errors  were  of  such  large 
magnitude  at  the  juncture  point  that  a properly  converged 
solution  could  not  be  obtained. 

Thus  far  it  has  been  shown  that  the  present  computa- 
tional approach  yields  good  numerical  solutions  in 
comparison  with  the  inviscid  theory  [8] . Further  comparisons 
are  now  presented  with  the  available  experimental  data 
for  spherically  blunted  cones.  Figure  16  presents 
surface  pressure  for  7.5°  half  angle  spherically  blunted 
cone  at  free  stream  Mach  number,  M =13.41,  free  stream 
Reynolds  number,  Rea  = 1515,  wall  to  stagnation  temperature 
ratio,  Tw/Tq  = 0.0741  and  a free  stream  temperature, 

T^  = 200 °R  corresponding  to  the  experimental  data  of 
Pappas  and  Lee  [19].  Significant  differences  from  that  of 
the  inviscid  results  are  noticed  at  the  sphere/cone 
juncture  point  for  this  case.  The  inviscid  pressure 
distributions  shown  earlier  through  Figures  11-13  predict 
a discontinuity  in  the  pressure  gradient  at  the  sphere/cone 
juncture  point.  However  for  the  present  low  Reynolds 
number  of  1515,  viscous  effects  smooth  the  discontinuity 
completely.  The  result  of  the  present  calculations  are 
seen  to  compare  well  with  the  data  shown  in.  Figure  16. 
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Note  also  that  similar  numerical  calculations  which  do  not 
include  the  proper  jump  effects  at  the  juncture  point 
are  seen  to  yield  seriously  erroneous  results.  Figure  17 
presents  the  ratio  of  the  wall  to  the  stagnation  point 
heat  transfer  for  the  same  test  case.  The  comparison 
of  the  present  calculations  with  experimental  data  when 
the  jump  effects  are  included  is  seen  to  be  excellent. 
Again  large  errors  are  observed  when  these  jump  effects 
are  ignored.  Figures  16  and  17  also  show  the  numerical 
results  obtained  by  Miner  and  Lewis  [5]  for  the  same 
body  (with  an  artificially  smoothed  juncture  point) 
under  identical  test  conditions.  It  is  seen  from  these 
figures  that  the  present  results  tend  to  show  better 
agreement  with  the  experimental  data  when  shock  jump 
effects  are  directly  accounted  for  at  the  sphere/cone 
juncture  point. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 

An  analysis  of  the  physical  flow  behavior  at  the 
sphere/cone  tangency  point  has  been  made.  This  analysis 
indicated  that,  independent  of  the  choice  of  the 
coordinate  system,  inviscid  theory  would  predict  a 
discontinuity  in  the  flow  gradients  only  at  the  surface 
at  the  sphere/cone  juncture  point.  However,  following 
the  analysis  of  Messiter  [7] , it  was  found  that  in 
the  limiting  case  of  very  high  Reynolds  number  this 
discontinuity  would  be  smoothed  out  by  the  sublayer 
interaction  effect  within  the  inner  scale  length.  It 
has  also  been  shown  that  the  use  of  a surface  coordinate 
system  introduces  discontinuities  in  the  flow  gradients 
relative  to  surface  distances  everywhere  across  the 
shock  layer  and  at  the  body  surface  at  the  juncture 
point.  Within  the  viscous  layer  this  gradient 
discontinuity  due  to  the  coordinate  system  would  tend  to 
vanish  to  lead  order  as  the  Reynolds  number  tends  to 
infinity.  Analytical  jump  conditions  were  developed 
at  the  sphere/cone  juncture  point  for  these  discontinuous 
flow  gradients  associated  with  the  choice  of  surface 
coordinate  system.  Finite  difference  formulations  were 
then  developed  that  account  for  these  embedded  gradient 
discontinuities  in  order  to  eliminate  numerical  difficulties 
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in  the  solution  of  the  full  viscous  shock  layer  equations. 
Such  solutions  were  obtained  by  a numerical  scheme  which 
utilizes  a time  dependent  relaxation  technique  for  the 
bow  shock  wave  shape.  Comparisons  of  the  present  results 
with  inviscid  solutions  at  high  Reynolds  numbers  and 
experimental  data  at  intermediate  ones  were  found  to  be  good. 

While  the  present  technique  for  treating  the  sphere/ 
cone  juncture  region  has  been  shown  to  yield  good  results, 
the  present  numerical  scheme  which  is  essentially  a time 
dependent  relaxation  technique  encountered  certain  diffi- 
culties worth  noting.  One  difficulty  that  does  arise  is 
the  oscillatory  behavior  observed  in  the  iterative 
solution  of  the  shock  layer  equation  at  some  point  down- 
stream on  the  surface.  While  an  under-relaxation  scheme 
was  found  to  effectively  remove  this  problem  for  most 
flow  conditions  of  interest,  it  is  recommended  that  in 
future  studies  the  continuity  and  normal  momentum  equations 
be  coupled  during  the  iteration  process.  Another  diffi- 
culty that  was  found  for  this  relaxation  scheme  was  the 
initialization  process  used  for  the  bow  shock  shape. 

While  this  technique  enjoys  a greater  degree  of  flexibility 
as  compared  to  previous  techniques,  non-the-less  difficulties 
were  encountered  for  the  low  cone  angle  cases  of  the 
spherically  blunted  cone  studies.  Future  studies  should 
consider  use  of  the  inviscid  bow  shock  shape  as  the  natural 
initial  shape  for  this  relaxation  technique. 
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FIGURE  4 MODEL  SHOCK  SHAPE 


AEDC -TR-77-20 


FIGURE  5 MODEL  SHOCK  SLOPE 
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APPENDIX  A 

ANALYSIS  OF  THE  FLOW  VARIABLES  IN  THE  JUNCTURE  REGION 

The  purpose  of  the  present  section  is  to  study  the 
analytical  behavior  of  the  flow  field  properties  across  a 
sphere/cone  juncture  point.  In  particular  it  is  necessary 
to  determine  what/  if  any,  restrictions  the  conservation 
laws  place  on  the  flow  variables  across  a curvature  discon- 
tinuity. For  simplicity,  only  two  dimensional  flow  is  considered. 

The  analysis  naturally  begins  with  the  integral  form 
of  the  conservation  laws  since  they  alone  are  capable  of 
accommodating  discontinuities  in  the  flow  variables  if  they 
are  called  for.  Since  the  viscous  shock  layer  approach  em- 
ploys approximations  to  the  differential  form  of  the  govern- 
ing equations,  it  is  first  necessary  to  identify  the  equiva- 
lent approximations  in  the  integral  formulation  of  the  prob- 
lem. This  is  first  performed  for  a fluid  element  located 
away  from  the  juncture  point  with  the  results  subsequently 
generalized  to  the  juncture  point. 

To  do  this, first  consider  the  infinitessimal  element 
shown  in  the  sketch  below  as  being  located  at  some  point 
s,n  away  from  the  point  of  surface  curvature  discontinuity. 

Note  that  from  the  geometry. 


As3  = 

(1  + icn  - 

An) 

A s 

(Al) 

ll 

w 

< 

( 1 + icn  + 

<An) 

AS 

(A2) 
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and  the  unit  vectors  required  in  the  conservation  equations 
are  given  as 


(A3) 
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For  small,  these  yield. 


t = :3  * 1/2  * zn  A*  + 


L = e _ - 1/2  k e As  + 
s — s n 


n, 


e_  = 


n 


en  - 1/2  k eg  As  + 


e_  + 1/2  k e As  + 
n s 


Also 


n- 


e®4  = es 


-*■  -*■ 

e„  = e_ 

n4  n 


The  continuity  equation  in  integral  form  is  given  as, 
//  (pv.n)ds  = 0 

which  can  be  evaluated  on  the  element  shown  to  yield, 

(P2u2  - p1u1)2An  + [p4v4as4  - p3v3As3]  = 0 

so  that  with  application  of  equations  (Al)  and  (A2) 
this  gives 


(p2u2  - p1u1)2An  + [p4v4(l  + Kn  + <An) 
- P3v3(l+Kn-KAn) ] As  = 0 


(A4) 


(A5) 


(A6a) 


( A6b) 
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so  that  as  As  -*■  0 


P2U2  = plul  <A7) 

The  differential  form  of  this  equation  can  be  recovered 
using  Taylor  series  expansion  to  write,  for  example,  that 

P2  ~ P 2 ps  ’ * * (A8) 

P1  = p “ TT  ps  + **•  (A9) 

which,  when  used  in  equation  (A6b)  gives 

(pu)s  + [{1  + kh) pv] n = 0 (A10) 

The  most  important  point  to  note  here  is  that  for  this 
equation  no  terms  are  neglected  in  the  thin  shock  layer 
approach  and  hence,  the  equivalent  thin  layer  version  of 
the  mass  conservation  law  across  a general  line  is  given 
by  equation  (A7) . This  is  not  the  case  for  the  momentum 
equation  as  shown  below.  The  momentum  equation  in  integral 
form  is  given  as: 


//  [(pv-n)v  + pn]ds  - //  t ds  = 0 


(All) 


Evaluating  the  first  (inviscid)  integral  on  the  four 
sides  of  the  control  volume  shown  in  the  sketch  yields 

T1  = -t<Plul  + Pj.>%  + Plulvl«n,]24n  +t‘»2u2  +P2,Ss. 

+ c2u2v2Sn,I24n  - [l>3u3v3Ss,  +(p3v3  +P3>VUs3 
^ J 3 


+ [P4v4u4es  +(p4v4  +p4)e]As 

4 4 


(A12) 
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which  is  now  rewritten  using  equations  (Al)-(A4)  to  give 
Ix  = {2An  [(p2  + P2u221  " (-Pi  + plui  *+  1/2kAs  (p2u2v2  + 

■> 

,plulvl)]  + As  Ep4u4v4(1  + icn  + icAn)-  p3u3v3(l  + ten  - <An)  ] )eg 

2 2 

+ {2An  [p2u2v2  - p1u1v1  - 1/2  kAs(p2  + p2u2  + Px  + >1 

+ As  t(P4  + p 4V4^ ) (1  + icn  + icAn)  - (p3  +p3v32)  (1  + wi-icAn)  ] }en 

Similarly  evaluating  the  viscous  term  yields# 


(t,  3, 


n 


n. 


x,  e ) 2An  + ( t 4 e + 
1s  S1  *s  s4 


n 


n 


) 6s, 


e ) 2An  + (t,  e_ 

2 s0  3 n, 

s 2 n 3 


AS, 


(A13) 


In  these  relations  those  terms  that  would  be  dropped 
in  the  full  shock  layer  approach  have  been  identified  with  . 
a single  underscoring#  while  those  additional  terms  whose 
contributions  would  be  dropped  in  a thin  layer  approach  have 
been  given  a double  under scoring.  These  terms  were  identi- 
fied here  by  first  proceeding  to  the  differential  form  of 
the  governing  equations#  as  was  done  for  the  continuity 
above#  marking  the  shock  layer  approximations  there  and  then 
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tracking  these  backward  to  their  source  terms  in  equations 
(A12)  and  (A13) . 

The  most  important  result  to  emerge  from  these  equa- 
tions is  that  across  a general  line  (i.e.,  As+0) , the  flow 
variables  in  a shock  layer  are  governed  by  the  conditions 
that 


p2  + P 2^2 

= px  + 

(A14) 

plUlVl  = 

P2U2V2 

(A15) 

Combining  the  second  of  these  with  the  continuity 
equation,  (A7)  requires  that  in  a full  layer  approach  the 
normal  velocity  be  continuous  across  a line,  while  in  the 
thin  layer  model,  no  such  restriction  is  encountered. 

The  same  procedure  can  now  be  applied  to  the  integral 
form  of  the  energy  equation  in  viscous  flows  which  is  given 
as 

2 ^ ’ 

II  p (h+  y-)  (vn)ds  - II  t • v ds  + II  n*q  ds  = 0 (A16) 

s s 

Evaluating  the  first  of  these  terms  on  the  four  sides  of 
the  element,  as  before,  yields, 

I3  = 2An[P2u2(h2  + u2/2  + V2/2)-?1u1(h1+uJ/2+v^/2)] 

+ As[p4v4(h4  + u4/2  + v4/2)  (1  + icn  + icAn)-  p3v3 
(h3  + u3/24-v3/2)  (1  + <n  - KAn)] 
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Similarly  the  viscous-  terms  become 


//  t.v  ds  = (-^  ux  + tx  v1)2An  + (t4  u4+t4  v4) As4 


n s 


s _n 


(A18) 


+ (t2  u2  + t 2 v2)  2An  + (t3  v3  + t3  u3)AS3 


n 


n 


and 


//  n*q  ds  = -q-[2An  + qA  + q?  2An  - q3  AB3  (A19) 


where  again  those  terms  equivalent  to  the  shock  layer 
approximation  have  been  identified  with  a single  under- 
scoring and  those  terms  equivalent  to  the  present  thin 
shock  layer  concept  have  been  given  a double  underscoring. 

It  is  now  possible  to  evaluate  the  resulting  constraint 
on  the  flow  variables  across  a line  by  setting  As+o  to  obtain 


h0  + u~/2  + v2/2 


= hx  + u‘/2 


(A20) 


Evaluation  of  this  relation  in  combination  with  equation 


(A7) , (A14) , and  (A15)  verifies  that  whereas  all  variables 
p,  u,  V(  and  h must  meet  certain  constraints  across  a line 
in  the  full  shock  layer  approach,  the  present  thin  shock 
layer  concept  provides  no  constraint  on  the  normal  velocity  v. 
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Now  that  the  equivalent  thin  layer  terms  have  been 
identified  for  a region  of  continuous  surface  curvature, 
one  can  proceed  to  the  sphere/cone  juncture  point.  For 
this  study  the  element  now  straddles  the  juncture  point 
and  interest  centers  on  the  relation  between  the  variables 
over  faces  1 and  2 as  Asa  and  As^o.  With  this  in  mind  it 
is  clear  that  the  only  terms  of  concern  in  the  conservation 
laws  are  those  integrals  over  the  end  faces  1 and  2.  Thus, 
for  example,  only  the  first  term  of  the  continuity  equation 
(A6a)  need  be  considered  and  thus  one  can  write  immediately 
that  across  the  juncture  point 

P2U2  = plul  (A21* 
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For  momentum  conservation,  equations  (A12  .)  and  (A13) 
represent  the  appropriate  point  of  departure  except  that 
now  one  cannot  employ  equations  (A1)-(A2)  and  (A4) . None- 
the-less  it  is  still  clear  that  only  those  terms  with  An 
as  a coefficient  are  of  interest  here  and  that  all  others 
can  be  ignored:  With  this  in  mind  one  can  proceed  to  seek 

the  limit  form  of  the  shock  layer  model  as  As  0 while 
keeping  in  mind  the  fact  that  on  each  side  of  the  juncture 
line  the  shock  layer  approximations  still  hold.  With  this 
approach  the  results  are  identical  to  those  presented  in 
equations  (A14)  and  (A15)  for  momentum  and  in  equation 
(A20)  for  energy. 
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APPENDIX  B 

■ DERIVATION  OF  SHOCK  DERIVATIVES 


The  shock  derivatives  dugh/ds,  dTgli/ds,  dpgh/ds 
and  dvgh/ds  are  derived  in  this  appendix  for  use  in 
the  viscous  shock  layer  solution. 

In  the  spatial  coordinate  system  the  shock  angle, 
a,  is  written  as,  (Figure  1) 

a = tan  ^ (c[x — ^ (®1) 

where  R = yfi  + ngcos<|)  and  xs  = XB  “ ngsin<|>  (B2) 


Hence  the  derivative  dct/ds  is  evaluated  as. 


da  = 1 d2R  dxs 

2l  dxs  ds 
s 


(B3) 


Note  that 

dxs' 

SI“ 


cos$(l+Kn_)  - n_  sin<|> 

5 S 


(B4) 


and 


dns 

31“ 


*=  (1+Kn_)  tan(a-iji) 


s 


(B5) 


combining  (B4)  and  (B5)  yields 

s _ t-t  \ cos  a 

as-  - (l+*ns)  — 


(B6)  - 
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substituting  (B6)  in  (B3)  and  noting  that  dR/dx  = tana 

s 

yields  after  certain  manipulations, 


da 

ds 


(l+Kn) 

O 


3 

COS  a d R 
cos  Cot— ) dx2 


(B7) 


It  is  to  be  also  noted  that. 


tana 


dR/ds 

dx  /ds 
s 


(B8) 


Hence  the 


second  derivative. 


2 2 
d^R/dx^ 
s 


/ 


2 2 
dzR/dsz 

(dxs/ds) 2 


d2xs/ds2  dR/ds 
(dxg/ds) 3 


can  be  shown  to  be 
(B9) 


Substituting  for  dxg/ds  from  (B6)  and  then  evaluating 
(B7)  yields  after  proper  manipulations. 


2 2 

da  _ d R , cos  ( ct— cf> ) ..  _ dR  ricsin  (2a-2<ti ) ., 
ds  dg2  (1+Kng)  cosij)J  ds  lcos<|)  (l+icns) J 


(BIO) 


Note  now  that  from  Reference  [16]  equation  (A-6)  gives 


du 


sh 


ds 


K1  ds  K K2 


(Bll) 


where 


K,  = (1 


" cos  (2a+B ) - sinacos  (a+B ) (^i)  K, 

Y *sh  V 1 

(B12) 
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I 

where  is  given  by  (A4)  of  Reference  [16]  and  is  of 
the  form, 


d_  fTsh* 
ds  (psh} 


K 


da 

1 3s 


(B13) 


Also, 


T 

K0  — cosa  costa+8)  + ^ — sina  sin(a+(3) 
2 y Psh 


(B14) 


Likewise  for  other  flow  properties,  [Ref.  16] 


dpsh  = K d0 
ds  *3  3s 


(B15) 


where. 


^3 = Ty+ry 


dT 


sh 


3T 


= K 


da 
4 ds 


(B16) 


and 


K = — - — s-  sin2a  + — - — * — j 
q (Y+l)  (y+1>  M*  sin ''a 


cosa 


dVsh  _K  da_K  K 
ds  K5  3s  K6 


(B17) 


where  Kg  and  Kg  are  given  by  expression  {A12)  of  Reference 

tl6l  • 
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APPENDIX  C 

CHARACTERISTICS  OF  THE  SHOCK  LAYER  EQUATIONS 

The  characteristics  are  obtained  from  the  inviscid 
full  shock  layer  equations  (17-22)  and  the  corresponding 
"strip  conditions"  given  by 

*>  = S ds  + rc  an  (cl) 


dv  = |v  ds  + |V 
3s  3n 


dn 


(C2) 


*»  = H *=• + !£  dn 
dp  - If  ds  + IS  3,1 


(C3) 

(C4) 


The  terms  in  the  normal  momentum  equation  (19)  which  are 
not  included  in  the  thin  shock  layer  version  of  the  full 
shock  layer  equations  will  be  marked  here  by  introduction 
of  a multiplicative  factor,  a,  which  would  be  zero  for 
the  thin  shock  layer  equations. 

Using  Cramer's  rule  to  identify  the  derivatives 

gives 

Is  = |A|/|B|  <C5) 

where 


99 


AE  DC-TR-77-20 


dn 

ds 

0 

0 

0 

0 

0 

0 

0 

0 

dn 

ds 

0 

0 

0 

0 

0 

0 

0 

0 

dn 

ds 

0 

0 

0 

0 

0 

0 

0 

0 

dn 

ds 

0 

P 

(1+Kn) p 

0 

(1+Kn) v 

u 

0 

0 

(1+Kn) p v 

pu 

0 

0 

0 

0 

0 

1 

0 

0 

a (1+Kn) p v 

apU 

0 

0 

(1+Kn) 

0 

0 

0 

0 

0 

2 

- (1+kh) va 

-ua2 

(1+Kn) v 

u 

du 

ds 

0 

0 

0 

0 

0 

0 

dv 

0 

dn 

ds 

0 

0 

0 

0 

dp 

0 

0 

0 

dn 

ds 

0 

0 

dp 

0 

0 

0 

0 

0 

dn 

ds 

-KpV 

P 

(1+Kn) p 

0 

(1+Kn) v 

u 

0 

0 

-ICpUV 

pu 

0 

0 

0 

0 

0 

1 

2 

KpU 

0 

a (1+Kn)p  v 

apu 

0 

0 

(1+Kn) 

0 

0 

0 

0 

0 

- (1+Kn) va2 

-ua2 

(1+Kn) v 

u 

(C6  ) 


Setting  | B | = 0 , expanding  and  simplifying  yields 

[udn  - (1+icn)  vds]  { otdn3u  (u2-a2)  - odn2ds  u2(l+Kn)v 

- adn2ds  v (1+Kn)  (u2-a2)  - dnds2u  (1+icn)  2 (a2-av2) 

“ dsdn  v(l+<n)u2a  + ds2dn  v2 (1+Kn) 2 ua 

+ ap2v2 (l+xn) 2 ds2dn  u - (1+Kn) 3ds3v (-av2+a2) } = 0 

(C7 ) 
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Note  that 

udn  - (1+Kn)v  ds  = 0 

gives  the  equation  of  a streamline  in  this  coordinate  system. 
With  further  simplification,  the  final  expressions  for  the 
characteristics  can  be  written  as 


1,2 


-uv(l+<n)  _ (1+Kn) 
2 + a 

/l,  2^  2 2 . 

/-(u  +av  -a  ) 

a 

(1  - u2/a2) 


(C8) 


Note  that  when  a=l,  which  corresponds  to  the  full  shock 
layer  equations,  the  characteristics  are  inclined  at  a 
Mach  angle  in  the  supersonic  flow.  However,  as  a approaches 
zero,  corresponding  to  the  thin  shock  layer  version  of 
these  equations,  the  characteristics  in  the  flow  field 
tend  to  become  perpendicular  to  the  surface.  This  result 
can  be  verified  by  beginning  with  the  thin  layer  equations 
and  repeating  the  above  derivation.  For  the  thin  layer  case 
it  can  be  shown  that  the  characteristics  are  given  as 


- ds2p2a2  (1+Kn) 2 [dnu  - (1+Kn)vds]2  = 0 (C9) 


Note  that  this  equation  indicates  that  either 

ds2  = 0 CC10) 

or  [udn  - (l+ien)vds]2  = 0 (Cll) 
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Since  (Cll)  represents  the  equation  of  a streamline  in 
this  coordinate  system,  the  equation  for  the  characteristics 
are  given  by  (CIO) . Along  these  characteristics,  compatibility 
conditions  must  be  satisfied,  these  being  obtained  from, 

2 2 2 2 
ds  p(l+icn)  I-ua  dn  +(1+Kn)va  ds]  [-p«u  dn+p  (l+icn)  udu 

+ (1+Kn)dp  + p (1+Kn)  icuvds]  = 0 (C12) 

Note  that  this  equation  (C12)  is  satisfied  along  a char- 
acteristic line  ds=0  indicating  that  no  other  additional 
condition  need  be  satisfied.  It  is,  therefore,  seen  that 
the  inviscid  set  of  thin  shock  layer  equations  predict 
coincident  characteristics  normal  to  the  surface  of  the 
body.  Davis  [20]  discusses  the  characteristics  and  the 
nature  of  these  equations  (i.e.  whether  they  are  elliptic, 
parabolic  or  hyperbolic)  when  viscous  effects  are  included. 
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APPENDIX  D 

ADI  FORMULATION  OF  S-MOMENTUM  EQUATION 


The  s-momentum  equation  is  written  in  the  following 
form  in  the  surface  coordinate  system. 


pu  us 
Tl+icnT" 


+ pvun  + 


jcuvp 

(l+«n) 


+ 


P s 

(l+icn) 


9 9 "1  9 n 

= [e  /(1+ienj  (r+n  cos$)  [ (1+ien)  (r+n  cos$)jt] 

n 

(Dl) 

where 

t = v [un  - ku/ (l+icn)  ] 

Using  the  transformation  given  in  equations  (5a~h)  this 
becomes 
2- 

TT  + “1  T?  + ”2  U + “3  + “4  35  0 tD2) 

3 n 


where  a^. 
Note  that 
and  a?. 


a-,  and  are  given  by  equations  (8a-d) . 


u 


sh 


and 


appear  in  the  coefficients 


From  Appendix  (B)  u ^ and  pg^  are  written  as, 


K, 


da 

ds 


(D3) 


Psh 


K 


da 
3 ds 


(D4) 
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where  da/ds  is  given  by  (BIO)  of  Appendix  (B) . Sub- 
stituting da/ds  in  (D3)  and  (D4)  and  then  evaluating 
the  coefficients  a2  and  yields 


d R . dR  J 

°2  y2  dg2  y3  ds  + y4 

d2R  , dR  . 

03  " Ys  ^2  + Ye  d¥  + Y? 

where, 

y4  = +a<  k2  + B + C + D 


(D5) 

(D6) 


(D7) 


y3 


+A  K Ksin(2a-2^) 

1 cos<|>  (1+Kn  ) 
s 


Y,  = -A  K 


cos  (ct-j>) 

1 (1+Kn  ) co s 4 

s T 


and 


A = + 


psh  nsh 


Kn 


sh  pu 


e P 


sh 


1+Kn  . n - 
sh  p 


B = - 


p shvshnsh 


n 


2 

e v 


sh  pv 


sh 


1+Kn  , n 

sh  p 


Kn 


C = - 


sh  pn 


1+Kn  . n 

sh  p 


Kn 


D = - 


_u  cos*  n , Kn  , 

/ sh  + sh>  sh 

(1+Knshn)  r+nghcos<|)n;  (l+Knshn) 


sh 


(D8) 

(D9) 
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In  a similar  manner. 


n 


Y7  - A,  pc  - 


1 ^ n 


sh 

sh 


* PJ 


y6  A1  ps^  a3  cos<|>  (l+Kn^) 


K 


K sin(2a-2j)). 


Yr  = A 


cos  Cot — 4> ) 


1 psh  U+icnBh)coB* 


and 


A,  = - 


psh  nsh 
~2 

e vsh 


n 


sh 


(1+Knsh,')  V U 


sh 


hence  equation  (D2)  can  now  be  written  as 


32 3u  3u  , , d2R  , dR  , 

+ ai  air  +h2  jpr + *3  a?  + y4)u 


2 - 
, , d R , dR  . v , 3U 

(Ys  Y6  di  Y7)+  “4  3£  = 


(DIO) 

(Dll) 

(D12) 


Further  rearrangement  yields. 


2-  - 2 

3 u , 3u  . , - , » d R . , - . »dR 

—2  * “1  3W  + (*2U  + V 77  +tT3u  + 7«) 

3n  as 


6;  ds 


(y4u  + Y7}  + a4  If  = ° 


(D13) 


For  the  first  half  time  step  of  the  present  alternating 
direction  implicit  method,  this  last  expression  is 
written  as 
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First  Sweep:  t = tn  + At/2 


2-*  - 2 n 

3 u * 3u  * f3  Rn  . 

an2  1 3"  2 3s2 

* * 3u  ~ 

+ S4  + 65  W = ° 


8R.]+  e*  3R 


n 


3 1 


’3  3s 


(D14) 


For  the  second  half  time  step  equation  (D13)  is  written 
as 

»%-4- 1 A 4- 

Second  Sweep:  t = t + 


32Rn+1 
9 s2 


- e. 


3R 


n+l 


3 1 


+ 8. 


3 R 


:n+1  r32U  3U 

3S  3^2  1 3n 


+ B, 


3u  i * 
3T  + 64] 


= 0 


(D15) 


where, 

61  = al 

® 2 = Y2a  + y5 

B3  = Y3G  + y6 
e4  = Y4u  + y7 

65  = o4  CD16) 


However  note  that  equation  (D15)  for  R must  be  independent 
of  n indicating  that  the  coefficients  of  this  equation 
must  be  independent  of  n. 
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It  can  be  shown  by  proper  substitution  that, 

83/62  “ -2k  tan(a-$)  (D17) 


whereas  by  using  the  first  sweep  equation 


+ »i  ft  + 6s  If  + Wi 

an 


2^n 


,n 


3“R“  . 2R  2R  . „ 3R‘ 


n 


3s 


T~  + It 


fit”  + 2k  tan (a-*)  3g 


(D18) 


Substituting  all  of  this  the  final  sweep  equation  is 
given  as, 


a2pn+1  3Rn+1 

2k  tan(a-<|>) 

3s 


3S 


fit 


,n+l 


S2Rn 

Ss2 


+ 2k  tan(a-$)  (2R  - Rn)  = 


(D19) 


/ 
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APPENDIX  E 


DERIVATION  OF  FINITE  DIFFERENCE  EXPRESSIONS 
AT  A JUNCTURE  POINT 


Consider  a typical  mesh  system  where  a juncture 
point  occurs  immediately  ahead  of  point  3.  Point  ’a' 
is  taken  to  be  immediately  ahead  of  the  juncture  point. 


Juncture  Point 
al 


1 2 3 4 5 

Jump  conditions  associated  with  the  first  and  second 
derivatives  at  the  juncture  point  are  given  in  the 
form 


<^£>  = 
dsz  a 


K (— ) 
K1  ds  3 


t^-4)  + K (||)  + K3 

dg2  3 2 ds  3 3 


(El) 


(E2) 


where  K1#  K2  and  K3  are  known  and  are  given  by  equations 
(41b)  for  the  model  problem.  Using  Taylor's  series 
expansion. 


R2  = 


- AR_ 


4i 2 

T~ 


a 


(E3) 


i A2  ' * A3  * ' 

r4  = r3  + ar3  + — r3  + <r  r3  + ••• 


(E4) 
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Rearranging  this,  (E3)  can  be  rewritten  using  (El)  and 
(E2)  as, 

R2  = R3  - R3  (K;l  - 2-  k2)a  + j-  k3  + §-  R3 

A3  ' ' 1 

- f-  Ra  + . . . (E5 ) 

so  that  finally, 

, (r4-r2)  AKi  , 

R = + £ + O(A^)  (E6 ) 

AtKj-  ^ K2+l)  2(1^-  j K2+l) 

Using  (E4)  and  (E5)  and  simplifying  yields, 

,,  (R2”r3)  + (Ki  " I k2)  (r4-r3)  k3 

T~  ^K1  ~ 7 K2  +1^  2 ^K1  “ 2 K2+1^ 


A ,(K1"  7K2)R3 
1 1 * 


f V I 


-R 


« i 


(KX  " J K2+1) 


-1  + O(A^) 


(E7 ) 


Expression  (E7)  gives  a formally  first  order  accurate 
finite  difference  form  of  the  second  derivative  immediately 
behind  the  juncture  point.  Note  that  the  first  order 
error  in  (E7)  can  be  estimated  for  the  simple  model 
problem  through  differentiation  of  the  differential 
equation  on  the  two  sides  of  the  juncture  point.  From 
the  model  equation  (38)  differentiated  and  evaluating 
at  "a"  and  "3",  one  obtains 
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Ra  + “la  Ra  + °2  Ra  = 0 


R,  + a,  + a„  = 0 


1-3  3 T 2 3 


(E8) 

(E9) 


so  that 


(K1  - § K2)R3"  - l“l3<Kl  - 7 K2>  - “lalRj' 


R3lK2  “la  " “2lKl  - I K2>+  “2  Kl] 


The  error  term  in  (E7)  can  now  be  estimated  as, 


K3  “la 


<E10) 


M“la(Kl_  I K2>-“la1  ••  I ' 4 

3 1 ±S_  R . |R  [ { o,  K - « <K  - | K,> 

3(Kl"  7 K2+1)  12  2 12  2 


+ e^K-^A]  + k3  ala  a)  J 3(XX-  j K2+l)  (Ell) 


Note  that  when  (Ell)  is  substituted  in  (E7)  and  formally 

■ the  second  order  accurate  finite  difference  form  of  the 

first  derivative,  dR/ds,  is  utilized,  this  would  result 

in  a formally  second  order  accurate  finite  difference 

2 2 

form  of  the  second  derivative,  d R/ds  . 

The  above  difference  formulations  are  valid  only 
when  a mesh  point  of  the  finite  difference  scheme  coincides 
with  the  juncture  point  where  discontinuities  in  derivatives 
are  encountered.  These  formulations  need  to  be  rederived 
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when  the  juncture  point  lies  between  two  mesh  points  of 
the  finite  difference  solution  scheme.  This  is  achieved 
in  the  following  analysis. 

Consider  the  mesh  points  as  shown  where  the  juncture 
point  lies  at  a finite  distance  of  away  from  point  2. 
Points  a and  b are  immediately  ahead  and  behind  the 
juncture  point. 

Juncture  Point 


Let  n represent  a fraction  of  the  step  size  A.  Thus, 


n = 5A  1 < 5 < 0 


(E12) 


Expressions  are  now  sought  for  the  shock  derivatives 
at  points  "2"  and  "3"  with  an  embedded  jump  occurring 
from  "a"  to  "b".  Using  Taylor* s series  expansion. 


R1  ~ R2  AR2  + T~  R2  “ 6~  R2  + 

1 2 ' 1 A ^ 3 

\ - r2  + iER2  + 2-  5 r2  + r-  E 

2 

R3  = Rfa  + A(l-e)R^  + J-  (1-5) 2 R^' 

• « » i * * ' 

Ra  = R2  + A^R2  + “1  R2  + ••• 


Ra  , = R2  + A?  Rj  + . . . 


(E13) 


Rj' ' + ...  (E14) 


+ 


3 


• i i 

*b 


+ . . . 


(E15) 

(E16) 

(E17) 
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After ' considerable  manipulations  these  give, 


U+A)!^  + R3  + (1-5)2K3-R2(5+A+1) 


.2  r£+A  0 
A [—5—  + y 


+ £A  + 


(1-0 


+ 0(A) 


■] 


(E18) 


and 

R2  = {R3-2R1P2-R2(1-2P2)+  ^(1-52)K3>/[A(P1+P2))+0(A2) 

(El  9) 


where 

Px  - S + A 

£2  2 
P2  = |-  + CA  + (1-0/2 


In  a similar  fashion  one  obtains, 

2 

R3'  = [R2+R4P1-R3(1+P1)“  J-  C2K3J/[ (P2+P1/2) A2)+0 (A) 

(E20) 

and 

2 

R3  = [2R4P2+R3(1-2P2)-R2+  ?2K3]/[ (p1+2p2)A]+0(A2) 

(E21) 


where 

P1  = ^(Ki~  7 CK2)  + (l-6)] 

P2  = [£(1-0  (Kj-  ? 5K2>  + 7~  + IC1”5)23 
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APPENDIX  F 


ERROR  ANALYSIS  OF  THE  "SHOCK  JUMP"  MODEL  PROBLEM 


For  the  simple  "shock  jump"  model  problem  studied  in 
Section  IV  the  exact  solution  is  known  and  thus  a detailed 
study  of  the  accuracy  of  the  finite  difference  scheme 
can  be  undertaken. 

First,  in  order  to  establish  the  truncation  error 
of  the  finite  difference  expression  (44b)  at  the  juncture 
point,  the  exact  values  of  the  shock  shape,  R,  were  used 

in  this  expression  to  obtain  a numerical  estimate  to  the 

2 2 . , . 
second  derivative,  d R/ds  , at  the  point  of  discontinuity, 

i.e.,  at  S=0 . 9 . Figure  18  shows  that  this  derivative 

linearly  approaches  its  exact  value  in  A,  indicating  that 

the  finite  difference  form  (44b)  is  indeed  first  order 

accurate  at  the  point  of  discontinuity.  Note  also, 

from  this  figure  that  at  any  other  point  such  as  S=0.5 

where  no  discontinuity  of  any  kind  is  present  the  finite 

difference  expression  (44b)  is  seen  to  be  second  order 

accurate  as  anticipated  since  in  such  a case  the  jump 

constants  1^,  K2  and  K3  take  values  of  1,  0 and  0 

respectively.  These  results  clearly  indicate  that  the 

local  truncation  errors  are  of  first  order  at  a jump  and 

second  order  everywhere  else.  It  is,  therefore,  evident 

that  the  second  order  type  behavior  observed  earlier  can 
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only  be  explained  through  a study  of  the  overall  truncation 
error  of  the  numerical  scheme. 

In  order  to  better  demonstrate  this  concept  a simpler 
problem  for  which  an  analytical  assessment  of  the  error 
is  possible,  was  considered.  Thus  we  take  the  simple 
problem  given  as 


b2R 


0 


Subject  to  the  boundary  conditions 
R ( 0)  = 0 and  R(l)  = 1 


(FI) 


(F2) 


which  has  the  exact  solution  as, 
bs  -bs 

R = t6  b - V-l  CF3) 

e - e 

i 

Consideration  is  now  given  to  three  different  schemes  for 
numerically  solving  the  same  problem,  the  first  two  of 
which  establish  the  method  for  assessing  the  overall 
truncation  error  and  the  last  of  which  directly  addresses 
the  present  problem. 

2 2 

Case  (a)  - Consider*  the  case  where  d R/ds  is 
represented  by  a second  order  accurate  expression  in  the 
entire  solution  region.  In  such  a case  it  is  possible 
to  show  straightforwardly  that  the  difference  solution 
approaches  the  exact  solution  (F3)  quadratically  in  As. 
Since  this  is  the  basis  for  later  studies,  it  is  shown 
below  in  detail. 
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The  difference  version  of  equation  (FI)  is  given  as 


*i+i  - Ei<2  + bV)+  Ri-i  - 0 


(F4) 


Ro  = 0 and  Rjj  = 1 (F5) 

/ 

which  has  a solution  of  the  form 

• j 

Ri  = A S1  + B S2  tF6) 


where 


sl,2  = 1+ 


2 2 
b"V 


4 


2 7 r— — 

b V . ,2.2 
2 — + b A 


(F7 ) 


Applying  the  boundary  condition  yields  the  final  solution 


(F8) 


In  the  limit  of  zero  step  size.  A,  the  difference  solution 
(F8)  should  approach  the  exaqt  solution  (F3)  quadratically 
in  A . To  verify  this,  consider  now  the  limiting  process, 
as  A -*■  0.  First  write  that 


sl,2  = 1 ± hL  + 


bV  + b_V  + 

2 “ 8 •** 


and  note  that 


S1  = S1 


s/A 


(F9) 


(F10) 
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which  is  now  rewritten  as 


bV 


1 

S1  - e 


s/A  log_tl+bA)  s/A  log  (1+  — ) 


{1  - 


25  * 5 . . 3b4A4  25b6A6 

D A + . . . . + ^ + 


128 


3b3A3 

8 

s/A 


8 


128 


.} 


Expanding  loge(l+x)  for  small  x and  rewriting  the 
resulting  expression  in  terms  of  exponential  functions 
we  have , 


S/A 

S1 


bs  -b2As/2  b3A2s/2 


b2As/2 


= e e 


. . . x e 


b6A5s/24 


3b3A3 


x e 


. e 


s/A  loge(l- 


{1+3bV  (1.  SbVj-1 


8 


8 


-b4A3s/8 

e 


Note  that  the  first  order  contribution  in  A,  is  precisely 
cancelled  out  yielding. 


s/A 

S1 


bs 

e 


CW3a2 

[(i+  + 


) (1+. 


) ...] 


It  is  seen,  therefore,  that  s®^4  approaches  its  exact 
bs 

solution,  e , quadratically  in  A . 

Due  to  symmetry,  a similar  analysis  for  the  remaining 
terms  results  in  the  following  forms. 
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s*  = ebs[l  + .AA2  + BA3  + ] 

S2  « ebs [1  + AXA2  + BXA3  + ...] 

s»  = eb[l  + aA2  + bA3  + 3 

S;  = eb[l  + axA2  + bxA3  + 3 

Hence  the  difference  solution,  as  A -*■  0,  can  be  written 
as 

(ebs-ebs)+A2 (Aebs-A1ebs) + A3 ( . . . ) 

R.  *=  — p- — l a L*  v ■ —5 — ■ " (Fli 

1 (e  -eD)+  A/(aeD-a1e°)+  AJ{...) 


which  can  be  manipulated  to  the  form 


+ 0 (A2) 


Figure  19  shows  the  computational  result  verifying  this 
analytical  derivation  where  the  function  "R"  and  its  de- 
rivative,  d R/ds  , are  seen  to  approach  their  exact 
values  as  a straight  line  in  the  square  of  the  step  size, 
A2. 

Case  (b)  - Consider  now  a case  when  equation  (FI)  is 
solved  with  an  algorithm  that  is  first  order  accurate 
in  the  entire  solution  region.  To  do  this  the  source 
term  R of  difference  equation  is  written  at  the  midpoint 
between  two  mesh  points,  point  a,  using  the  average  value 
of  R^  and  R^^  to  approximate  Ra.  A centered  second  order 
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difference  is  used  to  represent  . The  resulting  first 

order  accurate  difference  scheme  is  given  as 
i+1* 


i it 

i-l<  i 


777 


3 Ri+i{2"b2A2>  ~Ri(4+b2A2)  +2Rj,_1  = 0 (F12) 


7777 


Using  the  same  boundary  conditions  as  before, 


Ro  = ° , 


= 1 


(F13 ) 


the  difference  equation  solution  is  found  to  be 


i . x 

S1  "'s2 
•R  — 

i N N 
S1  " S2 


(F14 ) 


where  now 


„ _ •,  , w.  , 3b2A2  L 17b3A3  , 

= 1 + bA  + — £ — + — — + • • • • 


2 2 3 3 

_ _ , u.  , 3d  A^  17bJAJ  , 

s2  _ 1 - bA  + —4 32~  + . . . . 


Following  the  same  procedure  of  manipulation  as  before, 
it  can  be  shown  that 


i . bs  -b2As/2  b3A2s/3 

s^  = le  e e ' 


) 


/ 3/4  b2As  -9/16  b4A3s  . 

S • • • / 


Note,  here,  that  the  first  order  term  in  A,  does  not 
cancel  out  as  in  the  previous  case  and  the  final  expression 
can  be  rewritten  as 
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S^  = ebS[(l+^A£+  ...)(1+  ...)  ...] 

Figure  20  shows  the  computational  result  verifying  that 

2 2 

function  "R"  and  its  derivative  d R/ds  approach  their 
exact  values  linearly  in  A . 

Case  (c)  - Consider  now  the  case  when  one  point  in 
the  finite  difference  scheme  has  a first  order  error  while 
all  other  mesh  points  are  formulated  in  a second  order 
accurate  sense.  This  case  is  then  similar  to  the  original 
model  problem  of  Section  IV  where  the  local  truncation 
error  was  second  order  at  all  points  except  one  where  a 
first  order  local  error  was  encountered.  Through  study 
of  the  present  problem  one  can  more  easily  see  how  the 
introduction  of  a first  order  error  at  a single  point 
does  not  cause  the  global  truncation  error  to  rise  to  a 
first  order  level. 

For  this  study  consider  the  following  mesh  system: 


h-1  h h+1 

• • • 


where  in  regions  (A)  and  (B) , the  governing  equation  (FI) , 
will  be  written  using  a second  order  accurate  central 
difference  expression  while  at  the  mesh  point  i = h, 
a first  order  accurate  version  of  the  governing  equation 
(FI)  will  be  used  by  again  evaluating  the  source  term  R 
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between  two  mesh  points.  Note  that  this  procedure  is 
analogous  to  the  earlier  model  problem  (Eq.  38)  where 
at  the  model  juncture  point  the  difference  equation  was 
first  order  accurate  due  to  the  jump  effect  while  at  all 
other  points  it  was  second  order  accurate.  The  advantage 
in  the  present  simpler  problem  is  that  following  the 
procedures  of  cases  (a)  and  (b)  a compact  analytical 
error  analysis  can  be  made. 

The  difference  equations  to  be  solved  are  given  as 

Ri+1  " Ri(2+b2A2)+  R±_1  =0  i = 1,  ...  h-1  (F15) 


Ri+l<2-b2A2)-  R±(4+b2A2)+  2Ri_1  =0  i = h (F16) 


and 


*i+l  " (2+k2^2 ) Ri  + Ri_i  =0  i = h+1,  ...  N-l 


subject  to  the  boundary  conditions 


(F17) 


i = 0 R„  = 0 
o 


i = N 


= 1 


(F18) 


A finite  difference  solution  of  (F15)  between  mesh  points 


i=l  to  i=h-l  yields  the  following  result 


i l 
S1  “ s2 

Ri  “ Rh  '^TTTh1 
S1  S2 


(F19) 


whereas  between  mesh  points  i=h+l  to  i=N-l  one  obtains 

(F20) 


ri  - a si  + E s2 
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where, 


*h 


Sr 


Rh  sl, 


N 


sh  sN 

s2  sl, 


A = -^TT)/<s2  - -TT» 


S1  S1 


Sn 


S. 


N 


-h  _N 


B = (1  - 


®h  S1 » / , _N  S 2 ®1 


■)  / (s~  - 


) 


S , 


and 


Sl  = ! + b4  +^1+  feV  + .... 

,2  = 1.bAtbV.bV,  .... 


Note  that  both  solutions  (P19)  and  (P20)  depend  on 
which  is  still  unknown.  The  difference  equation  (F16) 
for  the  "h"  point  relates  R^  to  R^+1  and  through  the 

expression 

(4+b2A2) 

are  known  in  terms  of  R^  from  equations 
(F19)  and  (F20) , this  equation  would,  after  proper  mani- 
pulation, yield  an  expression  for  R^  in  terms  of  known 
quantities  such  as  b and  step  size  A.  It  is  then  possible 
to  analytically  assess  the  error  term  in  this  expression 
in  the  limit  of  zero  step  size.  Note  also,  here,  that 
this  error  would  propagate  to  other  mesh  points  through 
the  solution  (F19)  ahead  of  this  point  and  through  the 
solution  (F20)  behind  the  mesh  points.  It  is,  therefore. 
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necessary  to  first  estimate  the  error  in  "R^"  before 
an  attempt  is  made  to  analyze  the  errors  at  any  other 
points  ahead  or  behind  it.  R^  can  be  rewritten  in 
the  following  form 


.Vl^"  2 )+  Vl. 

(1+  b2A4/4) 


Substituting  R^+l  from  (F20)  and  R^-l  from  (F19) 
and  after  some  manipulations  one  obtains  the  following 
expression: 


V1- 


2 2 


K, 


s,  - s. 


s,  - s 


N 

s,  - s. 


N 


] = <- 


N 

s,  - s. 


9 k2a2 

2 ) (l-  K±~) 


N 


where, 

K, 


, N h h N N h+1  . N h+1, 

(s2  s1  - s2  s 1 s2  s1  + si  s2  * 


and 


K2  = (sj  sf1  - S*  s^> 


Consider  now  the  term 

h 


I = 


2 2 
b'V 


K, 


s,  - s. 


N 


- s. 


N 


which  can  be  reformulated  as. 


N 

,22  S™  l~s i 

i - M-  i4  t-i- 

2 Sj  s2rsl 


N 


) - -4  ( 


l-s. 


s,  - s. 


-±_)  ] {_± 1] 

s2-sl  - s* 
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Noting  that. 


1 - s. 


s2-  sx 


3 3 

.l.bi  bJAJ 
7 + 3 32“ 


and 


1 - s. 


s2~  S1 


1 bA  b3A3 
7 + T"  + ~32~ 


Rh  * 


A — % Ei-  [i+ 


N 


N 


N' 


’1  “ S2 


/ 2 ,1  . bA 

{~k  (i  + r- 

S r\ 


3 3 
bJAJ 

“3T- 


.N 


V ®1  , 1 . bA  . 

• • e t-  j + r-  + 

s. 


b343 

32 


+ . 


)}(- 


s,  - s. 


N 

s,  - s . 


N' 


+ . . . . ] 


Expanding  this. 


*h  = 


h h h 
S1  ” s2  + S1 


*2  bV 


N 


N 


S1  " S2 


N 

s,  - s 


N 2 


N 

S2  /I  . 

“E  (2  + 3“ 

So 


bA 


b3A3 

32 


SN  33 

S1  , 1 bA  bV 

”E  ' 7 + T~  + ~J7~  + 


)>  (- 


s,  - s. 


s?  - 


_N' 


s.  - s. 

X 4 

N 

s,  - s. 


bV 


bV 


N 


*{....}  + . 


The  extraction  of  the  error  term  from  this  expression  as 
A approaches  zero  is  a process  involving  further  manipula- 
tions, however  if  it  is  first  noted  that 
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s“  = e n [1+  A A + B A + . . . .] 

Where  A and  B are  known  constants,  and  that  a similar 
statement  holds  for  the  terms  Sj  , and  s?J  , then 
it  is  possible  to  rewrite  as 

bs,  -bs, 

_ h h 0 

*h  “ l^E ==l 1 + > 

e - e 


It  is,  thus,  observed  that  the  first  order  local  trunca- 
tion error  at  point  "h"  contributes  to  higher  order 
global  errors  and  the  difference  solution  still  approaches 
the  exact  solution  quadratically  in  A.  Note  also  that 
one  can  assess  the  error  in  ff*  from  the  differential 
equation 


ds  h 


Cb2  R)h 


so  that  on  substituting  for  R^,  this  finite  difference 
solution  gives 


=b2 

ds^  h 


] + b20 (A2) 


Hence 


ds  h 


+ b2  0 (A2) 

ds*  exact 
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indicating  that  the  finite  difference  second  derivative 
solution  would  also  approach  the  exact  solution  quadrati- 
cally  in  A . It  is  now  possible  to  direct  attention  to 
the  analytical  difference  solutions  in  regions  (A)  and 
(B) , ahead  and  behind  of  the  present  mesh  point  "h" . From 
equation  (F19)  it  is  directly  seen  in  the  light  of  previous 
derivation  and  the  fact  that  has  a second  order  error 
associated  with  it,  that  the  difference  solution  in 
region  (A)  would  approach  the  exact  solution  as  second 
order  accurate  in  the  limit  of  zero  step  size,  A.  The 
same  result  is  also  true  for  the  region  (B) , since  this 
region  is  written  as  second  order  accurate.  Figure  F-2 
shows  the  computational  result  for  this  case  where  a 
first  order  accurate  difference  equation  was  used  at  a 
point  s = 0.5  while  all  other  mesh  points  were  written 
in  a second  order  accurate  sense.  The  function  R and 
its  derivative  d R/ds  are  seen  to  approach  their 
exact  value  as  a straight  line  in  the  square  of  the  step 
size,  A,  verifying  the  present  analytical  result.  There 
remains  now  only  the  question  as  to  when  would  this 
first  order  local  truncation  error  in  the  numerical 
scheme  produce  an  explicit  first  order  global  error. 

This  point  is  addressed  through  Figures  F-3  and  F-4.  It 
is  seen  from  Figure  F-4  that  when  a computational  scheme 
utilizes  the  first  order  difference  equation  in  a fixed 
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region  of  the  entire  solution  regime,  i.e.  between 

0.3  <_  s 0.7,  the  function  "R"  and  its  derivative, 

2 2 

d R/ds  approach  the  exact  solution  at  a point  s=0.5 

linearly  in  A indicating  that  the  overall  truncation 

error  is  of  first  order.  Figure  F-3  shows  a similar  study 

where  the  first  order  difference  equation  was  only  used 

on  a fixed  number  of  points,  i.e.  (N)n  c - 3 < (N) A _ < 

U . 5 — 0.5  — 

(N) © 5+3,  in  the  entire  solution  regime.  This  figure 
shows  that  in  this  case  the  function  "R"  and  its 
derivative,  d R/ds  approach  the  exact  solution  as 
second  order  accurate  scheme  in  the  step  size  indicating 
that  in  this  case  .the  overall  error  is  of  second  order. 

It  is,  thus,  clear  that  a local  error  of  order  A will 
not  sum  up  to  a global  error  of  order  A if  it  only 
occurs  at  a finite  number  of  points  as  the  mesh  is 
refined.  The  global  error  will  only  rise  up  to  first 
order  level  if  an  infinite  number  of  points  contribute 
a first  order  local  error  as  the  mesh,  is  refined. 
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APPENDIX  G 


NUMERICAL  EVALUATION  OF  THE  JUNCTURE  POINT  JUMP  CONDITIONS 


The  first  and  second  time  sweep  of  the  numerical  scheme 

for  the  viscous  shock  layer  solution  has  been  discussed 

in  Appendix  (D)  as  represented  by  equations  (D14)  and 

(D15) . Both  of  the  time  sweep  solutions  require  informa- 

2 o 

tion  about  the  jumps-  associated  with  the  terms,  3 R/3s* 
and  3R/3s  at  the  juncture  point.  The  jump  condition 
associated  with  the  first  derivative,  dR/ds,  is  obtained 
straightforwardly  from  geometric  considerations  as 


fdR. 

'SeT' 

QS  sphere 


<1+V  <ff > cone 


(Gl) 


2 2 

However  the  jump  condition  associated  with,  d R/ds  must 
be  obtained  from  the  momentum  equation  (D13) , which  may 
be  rewritten  in  the  form 


u + a.u  + v.u  + y-  + a.u 

2*  tan<«-*>  i§  +[-^3 1.n  4 «_t, 

dS  <V2  u + YS> 


= 0 


(G2) 


Note  that  the  last  term  in  (G2)  must  be  independent  of  n. 
Hence,  evaluating  (G2)  on  the  two  sides  of  the  juncture 
point  yields, 
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(R">  sphere  ' 2tan(“^>  <E’>  sphere  + <re)  sphere  = 0 

(G3) 


<r' + (Fe)  = 0 
cone  cone 


(G4) 


where , 


Fe  = 


[unn+alun  4y4u  + ^7  + a4ug] 

(YoU  + Ye) 


(G5) 


These  two  can  be  combined  to  give 


(R  ) = (R  ) + 2tan(a~4>)  (R  ) , 

sphere  cone  Y sphere 


+ (Fe)  - (Fe)  , = 0 

cone  sphere 


(G6 ) 


The  term  "Fe"  should  be  a constant  across  the  shock  layer 
since  the  associated  equation  (G2)  is  independent  of  the 
normal  coordinate,  n.  This  was  verified  numerically  at 
every  stage  of  the  calculation  procedure  for  the  viscous 
shock  layer  code.  It  was  also  found  that  there  usually 
was  a point  in  the  shock  layer  where  the  denominator 
(y2u  + Y5)  in  equation  (G6)  would  pass  through  zero.  It 
is,  therefore,  obvious  that  at  such  a point  the  term  "Fe" 
would  be  in  error  due  to  numerical  truncation  process  and 
care  should  be  exercised  to  avoid  any  such  region  while 
evaluating  these  jump  conditions.  For  the  present  viscous 
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shock  layer  code  the  jump  conditions  were  evaluated  at  the 
first  grid  point  away  from  the  wall. 

Another  difficulty  that  was  encountered  in  the  numerical 

i 

evaluation  of  the  shock  jump  conditions  stemed  from  the 
manner  in  which  the  s-momentum  equation  {D2,  Appendix  D) , 
was  solved  in  the  present  form  of  the  viscous  shock  layer 
code.  The  sketch  below  shows  a typical  finite  difference 
mesh  configuration  for  the  present  scheme. 


Due  to  the  nature  of  the  ADI  algorithm  as  applied  here, 
it  was  necessary  to  solve  the  star  time  sweep  equations 
for  the  flow  properties  u,  etc.  at  the  numbered  points 
"1",  "2'',  "3"  etc.,  while  the  final  time  sweep  equations 
were  solved  for  the  shock  shape  at  points  Cq,  C^,  Cj,  Cg 
etc.  The  jump  conditions  given  by  equations  (Gl)  and 
(G6)  were  to  be  applied  in  the  middle  of  the  second  sweep 
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mesh  from  points  "a"  to  "b"  of  the  sketch.  This  then 
requires  that  the  driving  function,  Fe , defined  in 
equation  (G5)  be  evaluated  at  points  "a"  and  "b".  To 
do  this  properly,  the  value  of  Fe  at  "a"  was  obtained 
by  extrapolating  its  values  at  Cq  and  to  "a"  while 

the  value  at  ”b"  was  obtained  through  extrapolation  of 
Fe's  values  at  points  C2  and  . The  results  thus 
achieved  were  found  to  be  consistant  throughout  the 
calculations. 
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APPENDIX  H 

COMPUTER  CODE  FOR  THE  FULL  VISCOUS  SHOCK  LAYER 
EQUATIONS  FOR  SPHERICALLY  BLUNTED  CONES 

The  following  computer  code,  written  in  Fortran  IV 
was  used  to  obtain  numerical  solution  of  the  full  viscous 
shock  layer  equations  for  hypersonic  flow  past  spherically 
blunted  cones.  The  input  quantities  are: 

Main  Program 

DT  Time  step  size. 

IE  Number  of  mesh  points  in  the  n -direction. 

IEND  Number  of  mesh  points  in  the  s-direction. 

REYIN  Free  stream  Reynolds  number,  Re^. 

RMAC  Free  stream  Mach  number,  M . 

BO  Wall  to  stagnation  temperature  ratio, 

TEMP  Free  stream  temperature,  Tn  in  degree  Rankine. 

GAM  Ratio  of  specific  heats,  y- 

SIGM  Prandtl  number,  o. 

XFACT  Convergence  criterion  for  solving  the  governing 

equations  by  iteration. 

THETA1  Sphere/cone  angle  for  which  solution  is  desired,  6. 

THETA  Sphere/cone  angle  whose  solution  is  used  as  an 
initial  guess  on  the  shock  shape. 
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DTHETA  Increment  parameter  which  controls  the  increment 
in  angle,  A6. 

W Under  relaxation  parameter  required  during  a 

profile  iteration  procedure. 

WW  Relaxation  parameter  used  for  two  consecutive 

final  sweeps. 

J n-point  across  the  shock  layer,  where  shock 

jump  conditions  are  evaluated,  also  convergence 
criteria  put  for  profile  iteration. 

NJNC  Number  of  mesh  points  between  juncture  and 

stagnation  point. 

DY  Normal  step  size.  An. 

NITER  Number  of  profile  iterations. 

NTIME  Number  of  time  cycles. 

Input  Parameters 

THIN  Positive  when  thin  shock  layer  equations  used. 

THINI  Negative  when  full  shock  layer  equations  used. 

RUMP  Positive  when  jump  conditions  are  included. 

SWF AC  Positive  when  wall  slip  included. 

SSFAC  Positive  when  shock  slip  included. 

AHALF  Positive  when  an  initial  guess  for  half  the 

longitudinal  step  size  of  input  guess  needed. 
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MINF  Free  stream  Mach  number. 


TW/TO  Wall  to  stagnation  temperature  ratio, 

EPS  Defined  as,,  [y  («„  /Cp)/pB  a ]x/  . 

REY (INF)  Free  stream  Reynolds  number,  Ree . 

S £,  surface  distance. 


X Axial  distance  measured  from  nose. 


RSH  Shock  distance  measured  from  axis. 

NSH  Shock  stand  off  distance  normal  to  body  surface. 

XSH  Shock  axial  distance  measured  from  body  nose 

point. 

R Normal  distance  to  the  body  surface  from  the  axis. 

NSHP  dn  /d£ . 

S 

USH  u-component  of  velocity  behind  the  shock. 

VSH  v-component  of  velocity  behind  the  shock. 

TSH  Temperature  behind  the  shock. 

RSH  Density  behind  the  shock. 

PSH  Pressure  behind  the  shock. 


USP  dugh/dg. 

VSP  dvsh/d5* 
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TSP 

dTsh/d5- 

■ 

RSP 

dpsh/d5- 

PSP 

dpsh/as. 

PWALL 

Pressure 

* * *2 

at  the  wall,  p /p  U . 

00  00 

PW/PO 

Pressure 

ratio  at  the  wall 

CF 

if  it  it 

Skin  friction  coefficient,  2t  /p  U 

V7  oo  oo 

STAN 

Stanton  number,  q /(H  -H  ). 

w o w 

HEAT 

Wall  heat 

transfer,  q^/Cp!^3). 
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List  Of 

Name 

DERIV 

GEOM 

SHVALS 

MANISH 

PUSHPA 

PEQSO 

BOUND 

B0UND1 


Subroutines 

Function 

Calculates  the  derivatives  of  the  initial  shock 
shape . 

Calculates  body  geometry  for  any  given  longi- 
tudinal location's. 

Calculates  properties  behind  the  shock. 

Utilizes  shock  jump  conditions  to  evaluate  new 
shock  shape , final  sweep. 

Evaluates  further  shock  quantities  using  new 
shock  shape. 

Solves  tridiagonal  difference  equation. 

Provides  initial  coefficient  for  derivative 
boundary  conditions . 

Provides  initial  coefficients  for  derivative 
boundary  condition  at  s=0  on  the  shock  for 
final  sweep. 
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Flow  Chart  for  Sphere  Cone  Program: 


142 


AEDC-TR-77-20 


FORTRAN 

OCOl 

0002 


0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 
0013 
00  14 
0015 
00  16 
0017 


0018 

0019 

0020 
0021 

0022 

0023 

00  24 

0025 


0026 

0027 

0028 

0029 

0030 

0031 
00-32 


IV  G LEVEL  21 


MAIN 


DATE  ■ 76296  05/03/31 


1 

1 

2 

3 

4 

5 

6 

7 

8 
9 
1 

t 

1 


t 

2 

3 


IMPLICIT  REAL*8  IA-H,  0- Z I 
COM  MO  N/MA INI/ 

Tl(  201  1.T2I 20  1 ). TCI  201 1 .U1 1 201  ) .02(2011 . UC ( 20 l > . 

VI ( 20 1 1 • V2 ( 20 1 1. VC (2011.P1I201). P2 ( 20 1 1.PCI2011* 

T1N( 201 1 . T2N( 201  I ,TCN( 201  I .U 1 N ( 20 1 ) .U2N( 20 1 ) .UCN ( 20 1 ) . 

T INN ( 20 1 1 • T 2NN (201 ),U1NN( 20 1 1 . U2NN( 20 1 I .PCN( 201 > . 

AA(201  1 ,8B(201  1 . RV I SC ( 20 1 CON ( 20 1 1, VI SC (201 1 . 

RCON ( 20 1 I » CPST ( 20 1 >.RNSH(201 ).RCSF(201  I. 

Rl(201 1 . R2 ( 20 1 1. RC (201 l.P2N(20 11 ,PO( 20 1 1 .PON (201 1 . 

PS( 20 1 1 , V2N(201) .VS< 201 l.COl (201 }.C02(20 1 1 .PE ( 201 1. 

P1N( 201 1 ,P21( 201 1 .P22NI201 I ,P22<  201 1 , PFAC( 20 l 1 ,P21N(201  1 
•VG(201 l.VGN( 2011. P33( 201 1 .P33N(201 l.VGS( 20 1 1 . 
VO( 201 1 . V0N(201 1 .XM( 20 1 1 . PI TO( 20 1 1 .UCNN(  201  I 
COMMON  /PEOS/  OS  . DN  (201  I.  IM.IE.AK201  I.  A2  (201). A3(  20  11  . A4  ( 20 1 1 . 

X N( 202 1 


COMMON  / INSH/ 


DIMENSION 
DIMENSION 
DTMENSI ON 
O IMENSI ON 
DIMENSION 
0 IMENSI ON 
O IMENSI ON 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
COMMON/PUSHY/ 
COMMON /OUT SH> 

» 


CONO  . GAM  . S . UPSH 

EPS  . RMAC  . TPSH  * VISCO 

P2G ( 20 1 1 

P 1 3N( 202 1 
PF  AMI 2 101 
C12(2101.C11(210) 

CNS2P ( 40 1 .CNS2PP ( 40 1 
V1GI210  1 
V2G( 210 1 
CNS  2(110) 

VCDI (201.5) 

VCD 1( 201.51 

YNSH(llO)  . YNSP( 110 1 • Y NSPP (110) 

DER I VI  . THMAX 


XNS 


PPS 

• 

RRS 

• 

TTS 

• 

UUS1 

• 

vvs 

PPSl 

t 

RRS1 

• 

TTS1 

♦ 

UUS2 

* 

VVS1 

PSP 

• 

RRS2 

• 

TSP 

• 

USP 

• 

VVS2 

PPS2 

• 

RSP 

• 

TTS2 

• 

uus 

t 

VSP 

COMMON/BASU/  XS3  . CONE 

COMMON/KINNI/  XNSH( 1 I0).XNSP( 110 I.XNSPP(llO) 

COMMON/MANIS/  AXSH( 110) .AXSP( 110 ) . AXSPPI 1 10  > 

COMMON/MANU/EE1 .FF1 . I END. IEND1 . AAA)  (110) . AAA2( 110) . AAA3( 110) 
1.AAA4I 1101 

COMMON/CON/  NJNC.NJ1 , RUMP 

C0MM0N/MAIN2/  CNS.ALP.CONP  , AKK1  . ALP3 . PHI 3 . VVM . A KK2 .AKK3.X1SP 
1 .PHI .CSF2.RS2.TKK1 ,TKK2. TKK3 

COMMO  N/MN2/  RE  1 A , RE1 B . RE2A .RE28.RE3A.RE38 . YNPP J. YNSP J . YNSH J, DT . 
1 RSH1 . I. J.XNSP J.AMM3I 202) ,RSH 

DATA  BLNK  /•  '/.BNO/'NO'/ 

C COMPUTER  CODE  FOR  SPHERICALLY  BLUNTEO  CONE  USING  FULL  SHOCK  LAYER 
C EQUATIONS  (PROGRAMMED  BY  8 . N. SR  IV ASTAVA ) 

READ( 5, 30)  OT. THETA. THETA1 .WW.XFACT 
REAO( 5. 31 > J . NJNC. IE.  I END 
READ (5. 32)  RM AC . GO . RE Y I N. TEMP . GAM .S IGM 

30  FORMAT! 5F10. 6 ) 

31  FORMAT (4110) 

32  F ORM AT ( 6F1 2. 4 ) 

WRITE ( 6 .33 ) DT,  THETA , THETA1 ,WW . XF ACT , GA M. SIGM. J . NJNC 
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0033 

0034 

0035 
CO  36 
0037 
0030 

0039 

0040 

0041 

0042 
00  43 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 
00  52 
0053 
00  54 
00  55 

0056 

0057 

0058 

0059 

0060 
00  61 
0062 
0063 
CO  64 

0065 

0066 
0 0 67 
0068 

0069 

0070 

0071 

0072 

0073 

0074 
00  75 

0076 

0077 

0078 

0079 

0080 
0081 
0082 

0083 

0084 


33  FORMAT! 1H1 ,45X, 10HINPUT  DATA,/ 

17F 12. 6,2110) 

DTHET  A=3. 5000 
NT I ME  1 = 0 
N J1=NJNC-1 
NJ2=NJ1-1 
NJNCt=NJNC+l 

THMAX=THETA*3. 141 59225D0/1 80.0DO 
SMAX  = 3* 141 592653589793  2D0/2. ODO— THMAX 
AN  J l = N J 1 

DS=SMAX/1 ANJ1-0. 50D0) 

RUMP= 1.000 
R SHI  1=0.000 
THSL= 1. 0000 
THINI=-l.OD0 
AFULL— — 1 • ODO 
ALSL  = 1. OOOODO 
C0NVER=-1 .ODO 
36  CONTINUE 

THI N=THI N I 
35  CONTINUE 

TIME=0. ODO 
NTIME=0 
SWFAC— — l.OODO 
SSFAC=- 1.0000 

CALL  DERI VtDS.  I END , I END1 .YNSH.VN SP , YNSPP ) 

THMAX=THETA*3. 1 41 5922500/180.000 
SMAX=3. 1415926535  8979  32 CO /2 .ODO  — THMAX 
ANJ1=NJ1 

DS=SM AX/< AN J1 -0.50  DO ) 

NJ2=NJ1-1 
NJNC1 =NJNC+1 
77  CONTINUE 

IF!TIME.GT.35.0DO)  DT=!OO.ODO 
THMAX=THETA*3.  1 4 1 59225D0/1 8 0. ODO 
S MAX=  3.  141 592  65358979 32 CO/2. ODO  — THMAX 
DS=SMAX/<  AN J1  -0.50 DO) 

778  CONTINUE 

NTI ME l=NT IME 1+1 
DO  777  N= 1,1  END 1 

CNS2(  N)=YNSH(  N) 

CNS2P(N)=YNSP1N) 

CNS2PP1 N)=YNSPP(N> 

777  CONTINUE 

DERIV  1=— l.ODO 
XNSO— XNSH< 1 ) 

NTI  ME  = NT  I MEA’  1 
IM=IE- 1 
XN ( 1 ) =0. DO 
DO  15  N=  1 , I E 

IF!XN!NI.LE. l.ODO)  DY=0.0350DO 
I F! XN<N).LE.O.  64 9 9 999 DO ) DY=0.0150D0 

I F( XN(N).LE.O. 0 49  99999D0  > DY =0.001 ODO 


05/03/31 


144 


AE  DC-TR-7  7-20 


FORTRAN 

0085 

0086 
00  87 
0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 

0099 

0100 
0101 
0102 

0103 

0104 

0105 
<Sl06 

0107 

0108 

0 109 

01  10 
0111 
0112 

0 113 

0114 

0115 

01  16 
01  17 

0118 

0119 

0120 
0 121 

0122 

0123 

0124 

0125 

0126 
0127 
0 1 28 
0 1 29 

0 130 
0131 

01  32 
0 133 

0 134 

01  35 


IV  G LEVEL  21  ' MAIN  DATE  = 76296  05/03/31 

0N(N1=DY 

15  XN(N*1  > = XN( N)*DN(N> 

NI  TER=0 
RSH=0.0D0 
RSH1 =RSH1 I 
UUS0=0.0 
UR  SH=0 .0 
UPSH=O.ODO 
TPSH=0.0D0 
VISCO=O.ODO 
CONO=O.ODO 

ASL=1 . 2304*  C 2.  O-THSLI/THSL 
BSL=1 . 1750*12. 0-TMSL1/THSL 
C SL=2  »3071D0*(2» 0 DO— ALSL) / ALSL 
XNS=XNSO 
XNS1*XNS 
DS2=DS/2.0D0 
CK=1.0D0 
CSF=0.0D0 
SI  F=  1 • OD 0 
RS=0 • 0 DO 
RS2=0. 000 
X8=0 . 0 
CD  F=0 . 0 
CDP=0. 0 
CDP1 =0 . 

CDP2=0.. 

CDF 1=0 • 

CDF2=0. 

CDP0=0. 

CDFD=0. 

CNS=  (XNS1  + XNSJ/2. 

PO  IP=C  (GAM+1 . 1 *RMAC*RMAC  /2. 1 ** ( GAM/ 1 GAM-1 . 1 )/ C GAM*RM AC*RMAC* 

1 (2. *GAM*RMAC*RMAC/( GAM* 1. 1- (GAM-1. ) /(GAM+1.)|**{1. /(GAM-1. Ill 

TW=B0 *(  1 . 000/ ( ( GAM— 1 . 000 ) *RMAC*RMAC 1 + 0. 50 DO  » 

TB=TW*( (GAM-1.000 1*RMAC*RMAC* TEMPI 
CONP= 198. 6/(  ( G AM— 1 . 1 *RM 4C*RMAC *TEMP ) 

VI SRA=( 1 .O+CONP) *( 1 .0/( ( GAM-1,  1 *RMAC  *RMAC ) 1**1. 5/( 1. /( ( GAM- 1 . »* 

1 RMAC4RMAC l+CONP) 

EPS  = 1.0  / DSQRT ( REY IN* VISRA ) 

CALL  SHVALSd  . ODO . 0. ODO . 1 . ODO , 0 . C DO . TT SO .VVSO.UUSO ,PPS0 .1 1 
TT5=TT  SO 
DO  100  N=1 « IE 

RNSH( N)=CNS/(  1 . *CK*CNS*XN(N)  I 
RCSF(N»=CNS/(1.+CK*CNS*XN(N) 1 
U1 (N)=XN(N1 
U2(N>=XN(N) 

UIN(N)=1.0 
U2N(N)  =1  .0 
U1 NN( N ) = 0*  0 
UC(N)=XN(N) 

UCN(N)=1.0 
VI (N 1 =XN ( NI 
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0136 

0 t 37 
0133 

01  39 
01  AO 
0141 
0 142 
C 1 43 

0 144 

0145 

0146 

01  47 
01  48 

0 1 49 

01  50 
01  51 
0152 
01  53 
0 1 54 

0155 

0156 

0 1 57 
Cl  58 

01  59 
0 1 60 
0 161 
0 1 62 
0 1 63 
0 1 64 

0 165 
0 166 
0167 

01  63 

0169 

0170 

0171 

0172 

0173 

0174 

0175 

0176 

0177 

0178 

0179 

0180 
0181 
01-82 

0183 

0184 
01  85 
0186 
01  87 
01  88 


IV  G LEVEL  21  MAIN  DATE  = 76296  05/03/31 

V2 tNl=XN(Nl 
VC fNl=XN(NI 

T1 <N1  = 1. 0-t 1. 0-XNt Nl 1* t 1 . O-TW/TTSOI 
T2(N)=Tl (N) 

T1 N{ N>  =1 .O-TW/TTSO 
T2N( N)=T  INtNl 
T 1 NN(N 1=0.0 
TC(N)=T1 ( N 1 
TCN(N)=T1N<N1 

VISCINI=ITTS+C0NP1 *TCtNl **1. 5/<TTS*TC{N) 4C0NP1 

RVISCt Nl =( TTS»TC (N) +3. 0*CONP 1/(2.0*TC(N> * ( TTS4TC < N ) ♦CONP) 1 *TCN(N) 
CONlN>=VISC(N) 

RC0N(N)=RV1SC<N> 

PI (N)=l. 0 

P2CNI=1.0 

PC ( N) = 1 . 0 

PS<N)=0.0 

PO <N)  = 1.0 

PO  N{  Nl  =0 .0 

R1  <N)  = P1 (Nl/Tl  (N) 

R2<N)=R1 <N) 

RC(N)=R1 (N) 

PF  AC  (N  1=1.0 
PCN( N) =0.0 
PI N(  N 1 =0 . 0 
100  P2N( N 1 =0 .0 
AMI  )=0.0 
88(1 1=0.0 

VI  SCO=  ( 1 .0+C0NP1  *TTS**  1.  5/ITTS-+CONP  1 
CONO= V I SCO/S I G M 
00  5000  1=1. I ENO 
CRNI= 1.000 

IFtI.L7.601  W=0.  8 

IF<I.LT.91  W=0. 6000 

VI=I 

S=tri- 1. 0D014DS 

CALL  GEOMt S.DS2.RS2.CK2 .CSF2.SIF 2, XB21 
PHI  = DARCOSt  CSF2 I 
PHI 2= PH  I 

IFtI.EQ.il  PH  1 1 = 3.  141592  65  358979  32D0 /2. 0 DO 
PHIS2=(PHI2-PHI 1 !/DS*2.0D0 
IFtI.EQ.il  PHIS=PHIS2 

IFfCONE.LT. 0.0001  PHIS2=-1.0D0 

IFtCONE.LT.O.ODOl  PHIS=-1.000 

I Ft CONE.GT.O.  000 1 PHIS2=O.ODO 

I Ft  CONE .GT . 0.  ODO I PHIS=0.000 

0X0S1 =t  AXSPt I M-AXSPt I+I  1 1/2.000 
X1SP= XNSPt I 1 

XNSPM=t  XNSPt l 1 ♦ XNSPt 1 + 1 | 1/2.000 

ALP=DAT  AN( f VNSPt  I 1+YNSPtH-l  1 1 / t 2. 0D0*DXDS 1 1 1 

IFtRUMP.LT.O.ODOl  GO  TO  305 

I Ft  1. EQ.NJ1 !XNSPM=  t3.0DO*XNSPlNJl 1-XNSP  tN J21 1/2.0  00 
I Ft  1. EQ.NJ1 1 AXSPJ=t3. ODO* AXSP t NJI 1-AXSP1NJ21  1/2.000 


146 


AEDC-TR -77-20 


FORTRAN 

0169 

0190 

0191 

0192 
01  93 

0194 

0195 

0196 

0197 
01  98 

0199 

0200 
0201 

C202 

0203 

0204 

0205 

0206 
0207 
0200 

0209 

0210 
0211 
0212 

0213 

0214 

0215 

0216 
0217 

0216 

0219 

0220 
0221 
0222 
0223 


0224 

0225 

0226 
0227 
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MAIN 
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I FC I.EQ.NJ1 ) 
IFU.EQ.NJ1  > 

I F ( I.EQ.NJ1 | 

IF( 1.E0.NJ1 ) 

I FC  1.EQ.NJ1  I 
305  CONTINUE 

I F< I.  EQ.l ) 

SP  = DSIN(ALP) 

CP  * DCOSIALP) 
SPB=SP*SIF2*CP*CSF2 
CPB=CP  *S IF  2-SP*CSF  2 
2000  CONTINUE 

CALL  SHVALS  ( SP . CP 


YNPP J=(3. 000* YNSPPI N J l 1— YNSPP1NJ21 1 /2.000 
YNSHJ=( 3*000  *YNSH(NJ1 J-YNSH(NJ2» 1/2.000 
XNSH J=(  3. OOO*  XNSH(N J1  1—  XNSH ( NJ2 1 1/2.000 
YNSP J— 1 3.0D0*YNSP( NJ1 1-YNSPINJ21 1/2. 000 
ALP=OATAN( YNSPJ/AXSPJ) 


ALP3 (22*000/14*  ODO*DA  TAN( YNSP ( I+1I/AXSP1  1*1 1 I 1/2.0 


SPB.  CPB.  TTSH • VRSH  * URSH,  PPSH."  21 


C******* ************* ********************** 

IFII.EO.il  ALP3=  3. 1 4 1 5926535897932D0 /2. 000 
IFd.EQ.ll  PHI3=  3.  1 41592653569793200/2. ODO 
BLP*ALP3 
PHP=PHI 3 
AK=TTS/PPS 

AR=2.  ODO*GAM*RMAC*RMAC*DSIN(BLP 1**2.000-1  GAM- 1. 0001 

AR1=2. 0D0*GAM*GAM*RMAC*RMAC*DSIN(2.0D0*BLPI/( AR*C  GAM  + 1 .000 )** 1. 0 
100  1 

AP2=4 • ODO  * GAM*  OCOS ( BLP I / I ( GAM* 1 . 000  I *RMAC**2. 000*AR*DS I NC BLP1 * * 


13. 000  I 

AR3  — 4 .000* GAM*  *3.000*RM AC** 4.  000*DS 1N( 2 .000 *BLP I *DSI N ( BLP 1 **2.000 
l/( AR**  2. 000* ( GAM* 1 . ODO  1 ** 1 .000 ) 

AR4=2.0D0*GAM*RMAC**2. 0C0*DSIN( 2.0D0*BLPI/( AR*ARI * ( GAM*{ GAM*1 . ODO 
1 1/ ( GAM— 1 .ODOl  — 2.0D0*GAM* (GAM-1 .000 l/IGAM  + l. OOO I 1 
AR  5=4.  0 DO* GAM* GAM«OS IN( 2.000  *BLP l/t (GAM* 1.0D0)*AR*AR*OSIN(BLP)**2. 


10001 

ARRlsARl *AR2— AR3— AR4+ AR5 

AKK1  = -OSIN(2.000*BLP-PHP1*( 1 .0D0-( GAM-1 . ODO I* AX/GAM  I 
1 +DSI N { BLP I *DS I N(  BLP-PHP 1 * ( GAM- 1 . OOO 1 /GAM*ARR 1 
AK 1 1 =* DC  CS ( 2 . OOO *BLP-PHP I •(  1 .0  00-( GAM-1 . ODO 1 *AK/GAM 1 
1— DSIN( BLPI *0C0S( BLP— PH P 1 * ( GAM— 1 . 0 DO  1 /GAM*ARR 1 
AK  K2*0C0S( BLP 1 *DS I N ( BL  P — PHPI — ( GAM— 1 . ODOl /GAM*AK*OS I N ( BLP 1 * DCOS ( BLP 


1- PHP1 

AKKK2=-0C0S(3LP l*DCOS< BLP-PHP 1 - (G AM- 1 . 0 CO  I / GAM* AK*DS I N( BL F I *DS  I N 
1 (BLP-PHP 1 

AK6=2 • ODO  *GAM/ ( ( G AM* 1 . 0 DO  1 **2 . 0 00 1 *0S l N ( 2 . OOO* BLP I 

AK7=4 ,ODO*DCOS(BLP l/( ( GAM* 1 • OOO I **2 . OOO *RMAC**4 . 0 DC*DS I N ( BLPI ** 3. 
10001 

AKK6=AK6*AK7 

AKK 3* 2.000/1  GAM  + 1 . ODO I * OS IN ( 2. OOO* BLP I 
IFII.EO.il  GO  TO  302 

OALOS=(OCOS( ALP 3-PHI  3 1**2. 0 00/(( 1 .0  00  *CK*CNS I *DCOS ( PHI  3 1 1 1 
1*( YNSPP(  1 1 — < ( RSH— VNSH<  I I 1/DT 1*2.000 1-YNSP( I I *CK*0S I N ( ALP3* 2. ODO 

2— 2 • ODO *P HI 3 I / ( ( 1 .ODO*CK*CNSI *OCOS(PHI3I> 

302  CONTINUE 

I F( I .EO.l I 9 AL  DS  = ( (XNSPPIll  — (( CNS-XNSH( 1 1 1 /DT 1 * 2 . ODO I /( 1 . ODO* 
ICNSI-1 .ODOl 

USP=AKK1*DALDS*AKK2*PHI S 
V5P=AK1 1*0ALDS*AKKK2*PH  IS 


147 


AEDC-TR-77-20 


FORTRAN  IV  G LEVEL  21 


MAIN 


DATE  = 76296  05/03/31 


0228 

0229 

0230 

0231 

0232 
C 233 
C 2 3 A 

0235 

0236 

0237 

0238 
0 2 39 

02  AC 

0241 

0242 

0243 

0244 
C245 

0246 

0247 

0248 

0249 


0250 

0251 


0 2 52 
0 2 53 

0254 

0255 

0256 

0257 

0258 
0 2 59 
0 2 60 

0261 

0262 


PSP=AKK3*DALDS 

TSP=AKK6*DALDS 

RSP=<  GAM/I  GAM -l.ODO ) )*< PSP* TTS- TSP4PPS > /<  TTS*TTS) 

IFd.EQ.il  PSP=0.0D0 
IFIl.EQ.il  TSP=O.CDO 

IFII.EQ.il  V SP=0 • 0 DO 

37  CONTINUE 

TLP=ALP 
THP=PHI 
TK=TTS2/PPS2 

TR-2. ODO*GAM*RMAC*RMAC*OSIN(TLP  I*  DSINI TLPI-I GAM-1 .ODO I 
TRI  = 2 • 0O0*GAM*GAM*RM AC  *RMAC*DSIN(2.000* TLPI /(TR*t GAM+1 .ODO 1**1. C 
IDO  I 

TR2=4.0D0*GAM*OC0S(TLP)/{lGAM+1.0D0)*RMAC**2.0DC*TR*DSINITLPI* 

IDS IN<TLP)*DSIN(TLP! ) 

TR3=4.0D0*GAM**3.0O0*RMAC**4.0D0*DS INI  2. ODO* TLP|* DSINI TLPI **2.0 DO 

l/I TR**2. ODO* I GAM  + 1 .ODO 1**1.000  1 

TRA  = 2.0DO*GAM*RMAC**2.0  CO *DS I N I 2 .ODO *TL P 1 / I TR*TR I * I GAM* I G AM+ 1 .000 
1 1/ < GAM  — 1 .OD0)-2.0D0*GAM* <GAM-l,0D0)/lGAM+l.OD0)l 
TR5=4.  ODO*GAM*GAM*DS INI  2.0DO*TLP I /I  I GAM+ 1 . ODO ) *TR* TR *DSI N I TLPI **2. 
10DC1 

TRR 1 = TR 1 + TR2— TR3— TR4  + TR  5 

TKK1=-DS INl2.0D0*TLP-THP  1*1 1 .0  DO- < GAM- 1 . ODO I +TK/GAMI 
1+DSINt  TL P 1 * D S I N<  TLP-THP I * I GA M- 1 . ODO I /G AM *TRR 1 

TKll=+DCOSI2.0D0*TLP-THP)*(1.0D0-<GAM-1.0D0l*TK/GAM) 

1-DSI N( TLPI*DCO St  TLP-THP) *IGAM-1 . ODO I/GAM*TRR1 

TKK2=DC0SI TLP)*DSINITLP-THPI— (GAM-1.0DO)/GAM*TK*DSINITLPI*DCOSITLP 
1-thp  I 

TKKK2  =— DCOSI TLP)*DCOSITLP-THP|-IGAM-1.0DOI/GAM*TK*DSINITLPI*DSIN 
1 1 TLP-THP 1 

DALDS 1 = 1 DC OS  I TLP—  THP I **2.OD0/( I 1.0D0+CK2*XNS)*DCOS<  THP|)  I 
1*<  I YNSPPI I I +YNSPPI  I + l I 1/2. ODO- I I RSHl-I I YNSHI  I l+YNSHI 1+ 1 ) ) /2.0D0)  I 
2/DT) *2 .ODO l”CK2*0SIN( TLP* 2. ODO— THP*2.0D0 I/I l 1. 00C+CK2+XNS) * 

2DCQSI  THP ) I *( ( YNSPI  I l+YNSPI  1+1  1 1/2.000) 

I FIRUMP.LT. 0.  ODO)  GO  TO  306 
I Ft  I . EQ.NJ1 ) 

1 D ALOS 1 = < DCOSI TLP-THP I * *2. ODO/ ( C 1 .ODO  + CK  *XNS )*DCOS I THP ) I I 
2*1 YNPP J -{  IRSH1-YNSHJ) 

3/DT)*2 .ODO I - CK 2 *D5 INI TLP *2 . 0D0-THP*2. ODO I / I I 1. ODO+  CK2 * XNS I * 

ADC OS I THP) I *YNSPJ 
306  CONTINUE 

VSP1=TK11*DALDS 1+TKKK2*PHIS2 
C ********************************* ********* 

IFII.EQ.1I  VVM=VVS 
IFII.GT.l)  VVM= 1.000 

VI SCO=t 1 .O+CONP) *TTS**1. S/ITTS+CONP) 

CONO=V ISCO/SIGM 

REFAC=RRS* VVM*CNS/ I EPS* EPS* V ISCD  I 

VI S2=t  TTS2+C0NP)*T2I 1 ) * * 1 . 5/ ( TTS2* T2 I 1 l + CONPI 

XKSL  = V IS2*RRS*VVM*0SQRT( IGAM- 1 .ODO )*TTS2*T2( 1 )/GAM) / 

1 (PPS2*P2U  )*REFACI 

00  200  N = 1 . IE 
CPSTI N 1=1.0 
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0263 

0264 

0265 

0266 

0267 

0268 
0269 


0270 

0271 

0272 


0273 

0274 

0275 

0276 

0277 

0278 

0279 

0280 
0281 

0282 

0283 

0284 

0285 

0286 
0 287 
0288 

0289 

0290 

0291 

0292 

0293 

0294 

0295 

0296 

0297 

0298 

0299 

0300 
030  1 
0 302 

0303 

0304 

0305 

0306 


IF  IS. GE. 0.0301 ) GO  TO  160 

PFAC<N)  = 4.0D0*  <P2(N)+IPP S2/PPS3 -2 . 0 DO ) *P0 ( N ) )/(UUS2*DS> 

1 -XNSPI 2)*XN(N)*P0N(N)/(2.0D0*UUS2»CNS) 

GO  TO  200 
160  CONTINUE 
200  CONTINUE 

C SOLVE  ENERGY  EQUATION 

DO  500  N=l .IE 

A1 <N)=REFAC*VlSCO*CPST(N)*<UUS*XNSP( I)*RNSH(N)*RC(N)*UC(N)*XN(N> 

1 /( WM*CNS)-RC(N)*VCIN> ) / ( C ONO* CONI N I )+RCON<N  )4CK*RNSH<N> 

2 +RCSFIN) 

A4(N)=-REFAC*VISC0*CPST<N1*UUS*RNSHCN|*RC(N)*UC{N>/C VVM4CONO* 

1 CONINM 

A2  (N ) = A 4 < N ) ATSP/TTS 

500  A3  INI  = REFACFPPS*  V I SCO*  ( RNSHI  N I *UUS  *UUS*UC  < N ) *PF ACt  N > VVM*  VC  INI* 

1 PCNIN) ) /{ TTS*RRS*VVM*CONO*CON<  N) »♦ UUS*UUS*V I SC ( N » *V I SCO* 

2 (UCN<NI-CK»RNSH<NI*UC(Nn**2  /(  TTS*CONO*CON(  N ) I 
GAMP=GAM+1 .000 

GA  MM=G AM— 1.000 
RMACQ=RM ACTRMAC 
EPSQ=EPS*EPS 
SPQ=SP*SP 

FOGQ=4 . 0 00/ ( GAMPFGAMP ) 

OE  N=R  MAC  Q*RM AC Q*SPO 
CS l=SP*XNS/< EPSQ*C0N01 

CS2*-(  (URSH-CP) **2«-F0GQ  + GAM*SP0  + 12.0  DO /G AMH-FOGQ *G AMM ) /RMACO 
1 — FOGQ/DEN ) * 0.  5D0*SP*XNS/<  EPSQ*CONO*TTS2  > 

IF  (SWFAC)  501.501.502 

502  CBl=-t ./(CSL*XKSL» 

CB2=TW/(  TTS2*CSL*XKSL I 

CALL  B OU  NO ( T1NN.T1N.T1.CB1 .CB2 • E 1 . FI .CRN!) 

GO  TO  503 

501  El =0.0 

FI =TW/TTS2 

503  CALL  PE Q SOI TINN.TIN.T1 , T 2NN. T2N.T2.El *F1  .CRN I ,CS1 . C S2 , SSFAC  . 1 . 0 . I ) 
TT  S2G=  TTS2 

IF  < SSFAC)  521  .521  .522 

522  TP  SH  = T 2NI IE) 

TTS2=T2( IE ) *TT  S2G 

IF  (S.GE.0.0001 » GO  TO  525 
TT  SI =T  TS2 

525  TTS=(TTS2*TTS1 1/2.0 
DO  524  N=1 . 1 E 

T2NNI  N ) = T2NNt  N I FTTS2G/TTS2 
T2N(N)=T2N(N)*TTS2G/TTS2 
524  T2(N»=T2(N)#TTS2G/TTS2 

VI  SCO*! I .O+CONP ) *TTS**1. S/ITTS+CONP) 

CONO* V I SCO/SIGM 

refac=rrs*vvm*cns/ceps*eps*visco I 

GO  TO  523 
521  TPSH=0. 

C SOLVE  S MOMENTUM  EQUATION 

523  XU25=U2( 15) 
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0307 

XU2S=U2<J> 

0308 

VI  S2  = < TTS2+C0NP) *T2(1 1**1. 5/1TTS24T2C 1 1+CONP I 

0309 

XKSL  = VIS2*RRS*VVM*DSQRT(<GAM-1 .000  >*TTS2*T2(  1 l/GAM)/ 
1 ‘ (PPS2*P2(1  I*REFAC1 

0310 

DO  540  N=1 ,IE 

C 3 1 1 

R2 (N)=P21N1/T2(NI 

0312 

IF  (S.GE. 0.0001 1 GO  TO  541 

0313 

R1  (N  ) = R2  ( N ) 

0314 

T! (N>=T2(N) 

0315 

T1N(N1=T2N(NI 

0316 

T1NN(N)=T2NN(N> 

03  17 

541 

CONT  I NUE 

0310 

TC(N1=(T1(NI+T2<NI 1/2. 

0319 

TCN(N)=(T1N(N)+T2N<N1 1/2. 

0320 

VI SC1N >=<TTS+CONP) *TC(N>  ** 1 . 5/ ( TTS *T C 1 N 1 +CONP 1 

0321 

RVISC1 N| =( TTS*  TC<  Nl  + 3. 0*C0NP )/<  2. 0*TC( N) *( TTS*TC(N l+CONP ) > * 

0322 

CON(Nl=VISCCN) 

0323 

RCON(Nl=RVISC(Nl 

0324 

540 

RC (N)=PC(N1/TC(NI 

0325 

DO  600  N=1.IE 

0326 

A1  (N)  = RE FAC* ( UU5*XNSP(  I 1 *RNSH( N) *RC<  Nl *UC( N 1 *XN ( N 1 /( VVM4CNS 
1 -RC(N1*VC(N1 1/VISC{N1+RVISC(N1+CK*RNSH(N}+RCSF(N1 

0327 

A2  < N 1 = -REF  AC* (USP*RNSH(N)*RC(N1*UC(N 1 / VV M+ CK »RNSH ( N 1 *RC(N1 

1 *VC<  Nl l/VI  SC (N1-CK4RNSH1N1 *RVt  SC(N 1 - ( C K4RNSH ( N 1 + 

2 RCSF1NI 1 *CK*RNSH< N 1 

0328 

A3  (N 1 = -REF AC*PPS»RNSH(N 1 *PF ACINI/ ( V I SC (NI*RRS*VVM1 

0329 

600 

A4<N}=-REFAC*UUS»RNSH( Nl *RC( N) *UC(NI /( VI SCIN14VVM1 

0330 

CS 1 = SP  *SPB*UUS2  * XN  S/( EP  S »EPS* V I SCO* UR SHI 
1 -CK24XNS/I 1 . *CK2*XNS1 

0331 

CS2=-SP»XNS4(  CP  + VVS2*CPB1/CEPS*EPS* VISCOAURSH1 

0332 

IF  (SWFAC1  601.601.602 

0333 

602 

CB  1=-1 . / I ASL*  XKSL 1— CK24XN5 

0334 

CB  2-0. 

0335 

CALL  BOUND (U1NN.U1 N.U1  . Cel  .C02 . E 1 . F 1 , CRN  I } 

0336 

GO  TO  603 

0337 

601 

El  =0.0 

0338 

FI  =0.0 

0339 

60  3 

CALL  PEQSO (U1NN.U1 N ,U1  ,U2NN . U2N , U2 . E 1 ,F1  , CRN  I , CS 1 , CS2 , SSF AC . 

0340 

UUS2G=UUS2 

0341 

IF  (SSFAC1  621.621,622 

0342 

622 

UPSH=U2N( IE1-CK2*XNS*U2(  IE  1 / < 1 . + CK2* XNS 1 

0343 

UUS2  = U2<  IE  1 *UUS2  G 

0344 

IF  IS. GE. 0.00011  GO  TO  625 

0345 

UUS1=-UUS2 

0346 

625 

UUS= ( U US2+UUS1 1/2,0 

0347 

DO  624  N=1,IE 

0348 

U2NN(  N } = U2NN<  N 1 *U US2G/UUS2 

0349 

U2  N( N 1 =U2N ( N i • UU  S2  G/UUS2 

0 3 50 

624 

U2 (N | =U2 ( N 1 *UU  S2  G/UUS2 

0351 

GO  TO  623 

0 3 52 

621 

UP  SH=0  . 

C 

SOLVE  MASS  CONSERVATION  EOUATION 

0353 

623 

CONTINUE 
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0354 

0355 

0356 

0357 
0 3 58 

0359 

0360 

0361 

0362 

0363 

0364 

0365 

0366 

0367 
0 3 68 
0 369 

0370 

0371 

0372 
0 373 

0374 

0375 

0376 

0377 

0378 

0379 

0380 

0381 

0382 

0383 

0384 

0385 

0386 
0 3 87 
0 3 88 

0389 

0390 

0391 

0392 

0393 

0394 

0395 

0396 

0397 

0398 

0399 

0400 

0401 

0402 
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DO  640  N=1  • I E 

IF  (S.GE. 0.0001  I GO  TO  641 

U1  (NJ=U2(N) 

U1N(N1=U2N(N> 

U1NN(N|=U2NN(N1 
641  UC(N1=(U1(NI*U2(NI 1/2. 

UCNN(  N)  ={  U1  NN{  N)*U2NN( N J 1/2.0  00 
640  UCN(N) =(U1N(N1TU2N(N1 >/2. 

DO  6055  N=1.IE 

P2G(N1  = P2(N) 

V1G(N }=VC(N> 

V2GCN)  = V2(N1 
6055  CONTINUE 
XV10=P2(J1 
XVS0  = T2  <J) 

DO  700  N=2.IE 

AA  (N  ) = AA  ( N—  1 ltDN(N-l 1*(R2CN-1)*U2(N-11*R2(N1*U2(N1  1/2. 

700  BB(N)=BB(N-1  1 TO N ( N- 1 ) * ( R2 ( N- 1 1 *U2 ( N- I 1 *X N( N- 1 1 

1 +R2( Nl *U2( N) * XN( N ) 1/2. 

IF  (S.GE. 0.0001 1 GO  TO  705 
A! A=8. *BB( IE1 *RRS2*UUS2*CSF2/DS-DS 
BI B=4, *A A( IE1*RRS2*UUS2*RS2/DS-DS 
Cl C=-DS 

ROT=BI B*BIB-AIA*CIC 
XNS  = (-BIB  ♦ DSQR  T { ROT  1 1 / A I A 
X1R-XNS 
DO  701  N=l,IE 
XNR-XNS 
Cl  2(N1=C02-(N) 

C02(N1  =XNR*(RS2*AA(M+XNR*CSF2*BB(N1  1*R9S2*UUS2 

VC (N 1 = 8. *XNS*(  AA ( N > *RS2 ♦ XNS*BB < Nl *CSF2 1 * RRS2*UUS2/ ( OS*OS* ( 1 .♦ 

1 XNS*XN(NI 1**2  *RC( Nil 

V2(NI=VC(N1 

701  CONTINUE 
GO  TO  71  1 

705  A I A=BB ( I El *CSF 2*RRS2*UUS2 

B I 8=  A A ( IEI *RS2*RRS2*UUS2/2. 

CIC=-C01  { I E ) ♦( RS ♦CNS*CSF 1 *( ( 1 . +CK*CNS 1 *R RS*VVS-XNSP ( I 1 *RRS*UUS1*DS 

RQT=B I B*  BIB— A I A*C IC 

XNS  = (-BIB  * DSORT ( ROT  I 1 / AIA 

CIC1=-C1  1(IE1+(RS2  + X1R*CSF2)*(  (I ,*CK2*X1 R 1 *RRS2*VV  S2-XNSPM*RRS2 
1*UUS21 *DS 

ROT 1 = B1 B*BI B— A I A*C I Cl 
XN1 S= (—81 B+OSQRT( ROT1 1 1 /AIA 
DO  710  N=1.1E 
XNR=XNS 

C02(NI=XNR*(RS2*AA(N)+XNR*CSF2*BB(N> )*RRS2*UUS2 
VI*(C02( Nl-COl  (N1 1/DS 

VC (N1 =-V I/( RRS*VVM*RC(N1 * ( 1. ♦CK*CNS*XN ( N 1 ) * ( RS+CNS *XN ( N 1 *CSF 1 1 ♦ 

1 XNSP( I |*XN(N1*UUS*UC(N)/(VVM*(1.*CK*CNS*XN(N) 1 1 

X1R=XN1S 

Cl 2(N1=X1R*(RS2*AA(NI+X1 R*CSF2*BB(N} |*RRS2*UUS2 
VM  = ( C12(N1— Cll  (N! 1/DS 
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0410 

0411 

0412 

0413 

0414 

0415 
04  16 
0417 
04  18 
04  19 

0420 

0421 

0422 

0423 

0424 

0425 

0426 

0427 

0428 

0429 

0430 

0431 
0 4 32 

0433 

0434 

0435 

0436 

0437 

0438 

0439 

0440 

0441 

0442 

0443 

0444 

0445 

0446 

0447 

0448 

0449 

0450 


G LEVEL  21  MAIN  DATE  * 76296  05/03/31 

V2(N) =-VM/(RRS2*VVM*R2<  Nl*t  1.+CK2*X1R*XNIN>  ) * < RS2  + X 1 R*XN t N ) 4CSF2 I 
2I4-XNSPM4XNIN) *UUS2*U2 INI/I VVM*I 1 . 4-CK2 *X 1 R *XN ( N 1 I ) 

145  CONTINUE 

710  CONTINUE 

711  IF  tS.GE. 0.0001 ) GO  TO  715 
XNS1=XNS 

715  CNS=< XNSl*XNS>/2,0 

RSH1=RS24XNS*CSF2 
IF(I.EQ.l)  RSH1I=RSH1 

XSH=XB-CNS*SIF 
RSH=RS4-CNS*CSF 

IFITHIN.GT.O.ODOl  GO  TO  9116 

6020  CONTINUE 

DO  9117  N=1.IM 

VC(NJ=( W*V1GIN 1+1 1 .-Wl»VCINI 1 
V2INl=t  W*V2G(N)4-tl.  — W)*V2(Nll 
9117  CONTINUE 
9116  CONTINUE 

IF  I THIN. GE. 0.0  I GO  TO  716 
VVS2G=VVS 
VVS2G=VVS2 
VPG=VSP1 

716  CONTINUE 

DO  712  N=1 .IE 

RNSHINI=CNS/11.+CK*CNS*XNIN1 1 
IF  IS. GE. 0.000 1 1 GO  TO  713 
VI (N)=VCIN) 

PCSF(NJ  = CNS/(  1.  «-CK*CNS*XN(  N)  ) 

GO  TO  714 

71 3 RCSFIN  >=CSF*CNS/ I RS4- CNS*XN  I N I *CSF I 

714  CONTINUE 
VSINI=IV2(N)-V1(N| >/DS 

IF  I TH I N. GE. 0.0)  GO  TO  717 
VGI  N|  =V2(  N) 

VGSIN  1=VSIN1 

IFINTIME.EO.il  GO  TO  718 

I F(  I . LE  *2  I VGSIN ) = I VCOI IN ,1*1 l-VCDI IN, I I I /DS 
IFII.EQ.ll  VGSIN)  = (VC01 IN, 1 + 1 I -VCD I IN. I ) l/DS 

718  CONTINUE 

IFII.EQ.ll  V01NI=VCIN1 

717  CONTINUE 
712  CONTINUE 

DO  720  N=2, IM 

IF  I THIN. GE. 0.01  GO  TO  720 

VO  N(N1=I DNIN-1 i*V0 ( N+ 1 1/DNIN1-DNINI4V0 ( N- 1 l/DNI N- 1 I I / ( DNI N I ♦ 

1 DNI  N-l  I !♦ IDNINI-ONIN-l  I I * VO ( N I /< ONI N I 4DNIN- 1 I I 

VGNl N I =1 DNIN-1  )4VGIN+1)/DNINI-DN(N)*VGIN-1  I /DNI N-l  1 l/IDNINIt 
1 DNIN-1  I I 4-  IDNINI-ONI N- 1 I I * VGI N I /( DNI N I *DN I N- 1 1 I 

7 20  V2NINI=IDNIN-1  I * V2 I N+ I I / ON IN  I — DN I Nl * V2 I N — 1 I /ON  I N-l  I I /I DNI N 14 
1 DNIN-1 I }+(DNIN)-DNtN— II 1 ♦ V2 I N 1 / f DNI N 1 *DNI N- 1 ) 1 

IF  I TH IN . GE. 0 .01  GO  TO  725 

VON!  IEI=V0I1E1*(DNI  1M-1  14-2,4  DN  I I M ) l/IONI  IM  »♦(  DNI  IM  ) 4- DNI  I M- 1111 
1 -VO  I IE-1  1 *( DNI  IM-  1 l+DNt  IM 1 1/IDNt  I M > *DN l I M- 1 1 I 


152 


AEDC-TR-77-20 


FORTRAN 


0451 


0 4 52 


0 4 53 


0 4 54 
C 455 


0455 


0457 

0458 
04  59 

0460 

0461 
0 4 62 
0463 


0 4 64 


0465 

0466 

0467 

0468 
0 4 69 

0470 

0471 

0472 

0473 

0474 

0475 

0476 

0477 

0478 

0479 

0480 

0481 

0482 

0483 

0484 

0485 

0486 
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2 4V0< IE-2  > 40N<  IM)/<DN<  IM-1  >*<DN( IM)+DN< IM-1 ) ) > 

VGNt IE  > = VG< IEI 4<DN< IM-1  )+2,*DN( IMl  1/tDNC  IM I 4<DN< IM )4DN< IM- 1 ) ) ) 

1 -VG( IE- 1 ) 4<DN< IM-1 )4DN< IM) )/(0N( IM)*DN( IM-1 ) ) 

2 4VG<  IE-2  1 *DN<  IM)/<DN<  IM-1  >4<0N<  IM)4DN<IM-1  1 > > 

VON(  1 )=-VO<l  >4(DN<2)42.4CN<1  ))/(DN(I)*(DN(l  ) 4DNI2)  ) 1 

1 4 V0<2J4<DN< 2 )+DNI 1) )/(DN< 1 )*DN<2) ( 

2 -VO(3)*ON(1 >/{ DN< 2 >4< ON < 1 )+DN<2)) ^ 

VGN<  1 )=- VG< 1 > 4 <DN<  2)42. 4 0N< 1 1 > /( DN( 1 l*(ON( 1 I 4DN< 21  I 1 

1 ♦ VG<2J4<DN<  2 >4DN<  1 ) )/(DN< 1 )*DN<2)  ) 

2 — VG<3  J 4 DN  < 1 )/< DN<2J4<DN<1  )4DN<2> J J 
725  CONTINUE 

V2N< IE  J = V2< IE J4<DN< IM-1 )42.4DN< IM) !/<0N<  IM )*<DN<  IM )4DN< IM-I  ) J J 

1 -V2< IE-1 )4<DN< IM-1 )4DN< IM) )/<DN< IM)*ON< IM-1 ) ) 

2 4 V2< IE-2 ) *DN<  IM)/<DN<  IM-1  |4<0N< I M ) 4 DN  < IM— 1 ) ) ) 

V2N< 1 )=-V2<l J4(DN<  2)4  2.4 CN< 1 ) ) / < DN< 1 ) 4 < ON < 1 >4DN<  2)  ) ) 

1 4 V2< 2) 4 <DN<2J4DN< 1) )/<DN<l >*0N<2) ) 

2 — V2< 3) *0N  < 1 )/<DN<  2 )*<DN< 1 )4DN<2) ) ) 

C SOLVE  N MOMENTUM  EOUATION 

P21N< IE) =RRS2*UUS244  2 4CK24XNS/<  PPS2* <1. 4CK24XNS ) I 

P21<IE)=1.0 

P2N<  I E J — P2 1 N<  I E ) 

P2  < I E ) — 1 • 

RC  < I E)  = l .0 

IF  < TH IN.GE. 0 • 0 ) GO  TO  750 
IF<  I • EQ . 1 ) 

1 P3  3N < I E J *— RRS2 *V VS2  G* VVS 2G*< < 1 .-UUS24XNSPM/<  VVS2G4 < 1 . +CK2 ♦ XNS ) ) ) 4 
l VGN<  IE) 4UUS2*XNS4  VPG/VVS2G/ < VVS2G4  < 1 • +CK2  4 XNS ) ) I/PPS2 

IFd.GT.l  ) 

1 P3  3N  < I E ) =— RRS2  * < < VVS2G— UU52*XNSPM/ < 1 . ODO  4CK2*XNS) ) 4VGN< IE) 

24UUS2*  XNS/<  1 . 0004C  K 24 XNS ) * VG S<  IE ) )/PPS2 
P33<  IE  J=0.0 

P2N< IE )=P2N<IE )4P33N< IE) 

750  CONTINUE 
PC<IE)=1 . 

PCN<  IE )s<PlN< IE) 4P2N< IE)  ) /2. 

IF  <5. LE. 0.0001 ) GO  TO  755 

PFAC< I E > =! — XNSP < I ) *PCN < I E )/CNS4PSP/PPS  J/UUS 
755  CONTINUE 
PE  < I E ) — 1 . 

PS  < I E I =0 • 

R2< I E ) =1 . 

IF  < S. GE. 0.000 1 ) GO  TO  800 

CALL  SHV ALS< 1 . ODO. 0.000. 1 . ODO . 0. ODO . TTSO . VVSO ,UUSC , PPSO . 1 ) 

PON<  I E ) = 0.  0 

IF  <THIN.GE.0.0 ) GO  TO  760 
C PONI IE )=VVSG< 1 )4V0N< IE J/PPSO 

PON<IE>=  VVS  4V0N< IE J/PPSO 
760  CONTINUE 

PO  <IE)  — 1.0 
PI <1 E ) S1 . 

PI  NI  IE)  = P2N<  IE) 

Rt  < I E ) = 1 . 

VI  <IE)=1 . 
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0*87 

0488 

0489 

0490 

0491 

0492 

0493 

0494 

0495 


0496 


0497 

0498 

0499 

0500 

0501 

0502 

0503 

0504 

0505 

0506 

0507 

0508 

0509 

0510 

0511 

0512 

0513 

0514 

0515 

0516 

0517 

0518 

0519 

0520 

0521 

0522 

0523 

0524 

0525 

0526 

0527 

0528 

0529 


21  MAIN  DATE  = 76296  05/03/31 

KON=I  M 

DO  810  N=1,IM 

P2 1N(K0N(=RRS2*UUS2**2  *CK  2*XNS*  R2!  KON  )*  U2  ( K ON  ) **2  /( PPS2* ( 1 . +CK2* 
XNS*  XN(KON)  ) 1 

P2 1 ( KON  J =P21 ( KON+ 1 ) -ON ( K CN )* ( P2 1 N< KON  + 1 ) +P21 N< KON)  >/2. 

P2NI  KON) =P2 1 NIKON) 

P22IKON (=0.000 
P2 (KON )=P21 (KONI 
IF  ITHIN.SE. 0.  01  GO  TO  805 
IF( I .EQ.l ) 

1P3 3N(KONI=-RRS2*VVS2G*VVS2G*((R2(KON)*VG(KON>-R2(KaN)*U2(KON)*UUS2 

1 *XNSPM*XN(KON)/(VVS2G*( 1 . +C K2*X NS* XN( KON ) ) J )*VGN(KON) 

2 +UUS2*XNS*R2(KON)*U2( KON)*( VGS(  KON ) + VPG* VG ( KON ) / V VS2G ) 

3 /( VVS2G*( 1 .+CK2*XNS*XNIKON) ) ) )/PPS2 
IFII.GT.il 

1P33N(K0N )=-RRS2  * { ( R2 ( KON 1 *VG ( K ON 1-R2 1 KON ) *U2<  K0N)*UUS2 

1 *XNSP  M*  XN(KON>/(1.0 DO  * ( 1 .+CK2*XNS*XN( KO  N ) ) ) ) *VGN(  KON) 

2 +UUS2*XNS*R2(  KON)*U2(  KON) *( VGSI  KON)  ( 

3 /( 1.00041 1 ,+CK2*XNS*XN( KON) 1 ) 1/PPS2 

P33(  KON)  =P33(K0N+1  ) -DN<  K ON  )*  ( P33  N ( KON+  1 ) +P33MK0NI  )/2. 
P2N(K0N)=P2N(K0N )+P33N(K0N> 

P2 (K0NI=P2 (KON) +P33 (KON) 

I F( THIN.GT.0.0D0)  GO  TO  9222 
P2(K0N)=( W*P2G ( KON ) + ( 1. -W)*P2(KON) ) 

9222  CONTINUE 
805  CONTINUE 

R2 (KON)=P2(KON )/T2(KON ) 

IF  ( S. GE. 0.000 1 ) GO  TO  801 
P0N(KON)=0.0 

IF  (THIN. GE. 0.0 ) GO  TO  807 

PON(KON)=  VVS  *P0 ( KON ) *V0( KON) *V0N( KON ) /( PPS07T2 (KON) ) 

807  CONTINUE 

PO (KON)=PO( KON+1 )-DN( KON>* (PON(KON+1 ) *PON(KON) )/2. 

PI (KON )=P2 ( KON ) 

PI K( K ON ) =P2N (KON I 
R 1 (KON )=R2(KON) 

VI  (KON  ) = V2  ( KON  ) 

801  PE ( KON )=P21(KON)+P22(KON) 

PC(KON)=(P1(KON)+P2(KON) )/2. 

RC (KON )=PC(KON )/TC (KON ) 

PC  N(K0N)=(P1N(K0N) +P2  N(KCN) )/2. 

PS(K0N)=(P2(K0N)-Pl(K0N) )/DS 
IF  (S.LE. 0.000 1 > GO  TO  810 

PFAC(KON)=(PS(KON)-XNSP( I ) *X N( KON ) *P CN ( K ON ) /CNS+PSP*PC ( KON > /PPS ) / 

1 UUS 

810  KON=KON-l 

NI  TER  = NI  TER+  1 

IF(NI TER.GT.  300 ) GO  TO  6000 

TFACT1=XV10-P2( J) 

T FACT  2=  XV50-T2( J ) 

T FACT =X U25—U2 ( J ) 

I F(OABS (TFACT ) — XFACT ) 822,82  1,821 

821  GO  TO  2000 
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05  59 
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056A 
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822 

CONTINUE 

IF(DABS(TFACT11-XFACT1 

82A  ,821,821 

9777 

CONTINUE 

82A 

CONTINUE 

IF(CABS(TFACT21-XFACT1 

820,821,821 

820 

I F ( CQNV  ER . GT .0.000 1 
CCNVER=1.0C0 

GO  TO  823 

GO  fO  2000 

823 

CONTINUE 

IF  IS. GE. 0.0001)  GO  TO  830 
TST=(l.+TW/TTS01/2.0 

VI  SA=  I TTSO +C0NP1 *TST**1. 5/( TTSO*TST+CONP1 
VI  S3 -VI  SCO* VI SA 
REY=1.0/(EPS*EPS*VISC01 

XKXK-  (GAM-1.  l/IGAMU.  1 *( TTSO/TW+1.  1 *B0/( 2.*E PS*E  PS*VI  S3) 
IFIAFULL.GT. 0.0001  GO  TO  202 

GO  TO  222 
202  CONTINUF 
CNT-BLNK 
CWS-BLNK 
CSS-BLNK 

IF (THIN  .EC.  -1.)  CNT-BNO 
I F ( SWF AC. EQ  • -1.1  CWS-BNO 
IFISSFAC.EQ.  -1.)  CSS-BNO 
WRITE  (6,9221  CNT , CWS,  CSS,  IE,  I END,  OS 
922  FORMAT(1HO,3XA2«17H  THIN  SHOCK  LA YE° , 3X , A2 , 10 H WALL  SLIP,3X,A2, 

11 1H  SHOCK  SLIP, 5X, 15HN0  STEPS  IN  N «,IA, 

216H  NO  STEPS  IN  S =,IA,5H  CS  =,F5.3) 

WRITE  ( 6 , 92A 1 RMAC , BO,  EPS,  REYIN.  REY 
92A  FORMAT (1  HO, 5XAHM INF, 7X5HTW/TO , 7X3HEPS , 7X8HRS Y< INF! ,6X6HREY( SI / 
13F12.A.2E13.31 
222  CONTINUE 

830  REFAC*RP  S*VVM*CNS/ (EPS*EPS*V ISCO 1 

CFCH=2.*UUS*RRS*VVM*VISC ( 1 1 ♦ (UCN ( 1 1 -CK*CNS*UC ( 11 1 /REFAC 
HEAT-  TTS*P.RS*VVM*(CONO*CCN(  1 1 *TCN  ( 1 1 / VI  SCO«-UUS*UUS*VI  SC  ( 1 1 *UC  ( 1 1 
1 *UCN(1)/TTS)/REFAC 
STAN -HEAT/ ( 0.5* 1.0/( ( GAM-l. 0) *RMAC*RMAC 1 -TW 1 
XNSP ( I l-IXNS-XNSll/DS 
00  8A0  N-l.lE 

XM(N1=DSQRT( (UUS*UUS*UC(N1*UC ( N> ♦ VVM*VVM*VC (N 1 *VC ( N 1 1/ 

1 ((GAM-1.  1*TTS*TC(N) 1 1 

PO2P01-1 .0 

IF  (XM(N) .LE. 1.01  GO  TO  8A5 

P02P01- ( (GAK+1. 01  ♦XM(Nl*XM(NI/(2.*( GAM-1. 1 *XM( N1 « XM ( N 1 1 I ** 

1 ( GAM/IGAM-l. 1 1 /(2.*GAM*XM(N1*XM(N»/ (GAM«-1. l-IGAM-l. I /(GAM 

2 +l.ll**(l./(GAM-l.ll 

8A5  PIT0(N>=P02P01*PCINl*PPS«,(l.«-(GAM-l.l*XM(Nl*XM(N)/2.  )**(GAM/ 

1 (GAM-1. 11/POIP 

IFII.EO.il  VC(N1-VC(N)*VVS 
IFII.EO.il  V2 ( N1 -VC  C N 1 

IFlI.LE.5l  VCO I ( N» I 1-VC ( N) 

IFCI.EO.ll  VC0KN,I»«VC(N1/VVS 

IF(1.E0.21  VC01(N,I  l-VCINI/VVS 
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0575 
C 5 76 
0577 
C578 
0579 
0 5 SO 

0581 

0582 

0583 
0 584 

0585 

0586 

0587 

0588 

0589 

0590 

0591 

0592 

0593 

0594 

0595 

0596 

0597 

0598 

0599 

0600 
0601 
0602 

0603 

0604 

0605 

0606 

0607 

0608 

0609 

0610 
0611 
0612 

0613 

0614 

0615 

0616 

0617 

0618 

0619 

0620 
0621 
0622 

0623 

0624 


IV  G LEVEL  21 


MAIN 


DATE  = 76296 


05/03/31 


U1 (N)=U2(N> 

VI (N)  = V2 (N) 

T1  (N)  = T2(N) 

R1  CN ) = R2 IN) 

PI (N)=P2(N) 

T1N(N)=T2N(N) 

T1NN(N)=T2NN<N> 

U1 N(  N ) =U2N INI 
U1NNIN )=U2NN(N) 

Cll (N)=C12(N) 

840  COI(N)=C02CN> 

PWALL=PPS*PC(  1 ) 

IF  (SWFAC)  843  « 843  » 844 

84  4 PWALL  = PWALL-BSL*EPS**2*DSQRT( I GAM- 1 . ODO > /( GAM4TTS4TC (1)1 >*VISCO 
1 4V ISC< 1 1 *TTS*TCN( 1 1/CNS 

843  CONTINUE 

IF  IS. LE. 0.000 1 » GO  TO  841 
C0P2  = 4 . ♦ RS*S I F *P  W ALL 
C0F2=2.*RS*CSF»CFCH 
CDPD=CDPO»-(  CDP1+C0P2  )*DS/2. 

CDF0=C0FD»(CDF1 +CDF2  >ADS/2. 

CDP=CDPD/(  RS*RS ) 

CDF=COFD/( RS4RS 1 

841  IF  I S. GE.O. 000 1 ) GO  TO  842 
PW  ALO=PW  ALL 

CDF=0. 0 
CDP=2. 04PWALL 

842  COTOT=  CDF+CDP 
COPl=COP2 
CDF1 =COF  2 
PWRAT=PW ALL/PWALO 
XNS1=XNS 
UUS1=UUS2 
VVS1=VVS2 
TTS1=TTS2 
PPS1=PPS2 
RRS1=RRS2 


I L= I END  — 1 

IF( I. EQ.IENDl  XNS2 1 =RSH1 

I F( I . EQ . I END ) XNS20=RSH 

IF( AFULL.GT.O .000 1 GO  TO  203 

GO  TO  223 
203  CONTINUE 

WRITE  (6.926)  S . X8 . RS. CNS . XNSPI I ) . XS H. RSH. N I TER 
WRITE  (6.928)  UUS,  VVS,  TTS.  RRS . PPS 
928  FORMAT { 1 HO . 8X3HUSH  * 10  X3HVSH . 10X3HT5H. 10  X3HRSH. 10X3HPSH/3X6F13.6) 
WRITE (6.864)  USP.VSP.TSP.RSP .PSP 

864  FORMAT ( 1 HO , 8X3HUSP, 1 0X3HV SP . 1 0X3HTSP . 1 0 X3HRSP . 1 0 X3HPSP/3X6F 13.6) 
WRITE! 6.927 ) CFCH, HEAT, STAN, COF, COP. CDTOT, P WALL .PWRAT 
927  FORMAT ( 1H0.5X. 2HCF . 10X4HHE AT .8  X4HST4N ,8X  3HCOF .9X3HCDP , 9X5HCDT0T , 7X 
1 5HPW ALL . 7X5HPW/P0/8F1 2.6  ) 

926  FORMAT (1  HO ,SX, 1HS. 1 1X1HX . 1 1 X 1 HR , 1 1 X3HNSH , 9X4HNSHP , 8X3HXSH , 9 X3HRSH, 
I5X.7HN0  ITER/7F12. 6,16) 


156 


AEDC-TR-77-20 


FORTRAN  IV  G LEVEL  21 


MAIN 


DATE  = 76296  05/03/31 


0625 

0626 

0627 

0628 

0629 

0630 

0631 
, 0632 

0633 

0634 

0635 

0636 

0637 

0638 

0639 
06A0 

0641 

0642 

0643 

0644 

0645 

0646 

0647 

0648 

0649 

0650 

0651 
0 652 
0653 
0 6 54 

0655 

0656 

0657 

0658 
0 659 
0660 
0661 
0662 

0663 

0664 

0665 

0666 
0 6 67 
0668 

0669 

0670 

0671 

0672 
0 673 
0674 


556  FORMAT! 3X.3F10.6I 
GO  TO  223 

930  FORMAT (l HO,  BX5HU/USH . 8X 5HV/ VSH . 8X5HT/TS H. 8X5HR/RSH, 4X 1 1 HP/PSH1 APP 
1R1 * 6X5 HP /PSH , 8X5HN/NSH , 8 X4HMACH • 9X4HPI TO  1 
932  FORMAT (3X.9F13. 61 
223  CONTINUE 

IF1I.EQ.NJ11  AM=  1.  0D04XNS 
AAK3= 1 « OOO 

IF! I • EQ.NJ1 | A AK5=2» ODO*AM*DTANl ALP-PHI  1 

IF(I.EQ.l)  GO  TO  999 

A A3  = - YNSPPI 1 1 ♦2.000*CK*DTAN( ALP3-PHI 3 » * YNSP I I ) *2 . 000* ( 2.0 D0*RSH- 
1 YNSHI  I 1 l/DT 

AAAI C I )=-2,ODO*CK*OTAN< ALP3-PHI31 
A A A 2(  l l=-2.  OOO/OT 
AAA  3(1  1 = AA3 
A A A4 ( 11=0.000 

IF( I. EQ.NJ1 ) ALP3A=( 3. OOO* AAA3INJ1 l-AA A3INJI-1 ) ) /2.0O0 
IF(  I.EQ.NJNC1  1 ALP3B=(  3. ODO* AA A31 N JNCI - AA A 31 N JNC+ 1 1 1/2.000 
999  CONTINUE 
RS=RS2 
S-S+DS2 

CALL  GE0M(S,DS2.RS.CK,CSF,SIF,X8  1 
PHI  3=0 ARCOS(  C SF 1 
PHI 1 * PH 1 3 

1F( C ONE. GT. 0 . 0U0 1 PHiS=0.000 

IF(CONE.GT. 0.0001  PHIS2=0.000 

I F1CONE.LT. 0.000 1 PHIS=— 1 .000 

IFICONE.LT. 0.0001  PHI S2=— 1 . ODO 

ALP3= DATAN  1 VNSP (I+l 1/AXSP1I+1 11 
SP3=DSI N( ALP3 1 
CP3=OCOS( ALP31 
SP03=SP3*SIF+CP3*CSF 
CPQ  3=CP3*  SIF  — SP3*  CSF 
PHIS= (PHI3-PHI  I /OS  *2. ODO 
5010  FORMAT 16F10.4 1 
RS2=RS 
N1TER=0 

CONVER=— 1 . ODO 
IF1I.EQ.1I  CNSO=CNS 

5000  CONTINUE 

CN* 1 TIME  SWEEP  STARTS 

ENDNSH=(  XNS2  1 — XNS20 1 4-  XNS21 
CALL  BCUND1 1EE1.FF11 
EE1=0.000 
FF 1=0. OOO 
A AK4  = ALP3  B— ALP3 A 
I ENO 1 =1 END* 1 

CALL  MAN I SHI YNSPP . YNSP . YNSH . ENDNSH ,DS . AAK3.AAK4, AAK5.AM1 
CALL  PUSHPAl Y NSH.DS , IENO.YNSPI 
T IME-TIMETDT 

IF! AFULL.GT.O. ODO. AND. THIN. GT. 0.0001  GO  TO  6005 
IF  1 T HETA. GT . THET  All  GO  TO  8010 

IF! AFULL.GT.O. ODO 1 GO  TO  6005 
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0676 

0677 

0678 

0679 

0680 
0681 
0682 

0683 

0684 

0685 

0686 

0687 

0688 

0689 

0690 

0691 

0692 

0693 

0694 

0695 

0696 

0697 

0698 

0699 

0700 

0701 

0702 

0703 

0704 

0705 

0706 

0707 

0708 

0709 

0710 
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DO  8009  N=1,!EN01 

CONV=DABS<YNSH<N»-CNS2 IN)  ) 

I FICONV.GT.O.  001001  GO  TO  8010 

8009  CONTINUE 

CONV2=DABSCCNSO-XNSH< 1 I > 

GO  TO  18 

8010  CONTINUE 

00  78  N=1,IEND1 

YNSHI N)=WW*CNS2(N)+YNSHCN)*( 1.0-WW) 

YNSP( N}=WW*CNS2P( N I *YNSP(N>*( 1. 0-WWI 
YNSPP<N1=WWYCNS2PP(N| +YNSPP<N)*<1 .O-WWI 
78  CONTINUE 

CALL  PUSHPAIYNSH.DS, IEND.YNSP) 

IFINTIMEl .GT.2.AN0.THETA.GT.THETA1 1 THETA=THETA  -DTHETA 

I F<  NT IME1 . GT.2. AND.THET A. GT . THETA l ) NT I ME  1 = 0 

GO  TO  77 

18  CONTINUE 

IFIAFULL.LT.O.OOO)  NTI ME=0 

AFULL=1 .0D0 
GO  TO  77 

6005  CONTINUE 

IF(THIN.EQ.-1 .000 ) GO  TO  19 

1 END=  I E NO— 1 

TNI N=  TH  I NI 
NTI ME=0 
AFULL =-1.000 
GO  TO  77 

5050  CONTINUE 
19  CONTINUE 

00  1200  N— 1 . I ENO 1 

WRITE (7 ,1201 ) YNSH1N 1. YNSP1N) , YNSPP<  NI  » AXSP( N ) ,XNSH( N> . XNSP(  N) 
1200  CONTINUE 
1201  FORMAT! 6F12. 8 > 

6000  CONTINUE 
STOP 
ENO 
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0005 

0006 

0007 

0008 

0009 

0010 

0011 
00  12 

00  13 
00  14 
0015 
00  16 


00  17 

00  18 

0019 

0020 
0921 
0022 

0023 

0024 

0025 

0026 

0027 

0028 
0029 

00  30 

0031 


0032 


0033 


0034 


0035 

0036 


IV  G LEVEL  21  PEQSO  DATE  * 76296  05/03/31 

SUBROUTINE  PEQS0IW1NN.W1N.W1 , W2NN, W2N. W2 .E 1 .FI .CRN  I .CS1 .CS2.SSFAC 
1 .ENO.I  BACK) 

IMPLICIT  RE AL  4 8 (A-H.  0-Z) 

COMMON  /PEOS/  DS. ONI 201 ) • IM. IE. 41(201 I.A2I201 )• A3 ( 20 1).A4(201). 

1 X N ( 202) 

01  MENS  ION  MINN  I 201  ) .MINI  201) ,W II 20 1 > 

01  PENSION  EI201  ) . F ( 20 1 ) , W2NNI201  ) . W2NI 20 1 ) . W 21 20 1 ) 

61  1)=E1 
FI  1>  = F1 
DO  10  N=2,IM 

A=(2.0D0-A1( Nl 4DNIN) )/(DN(N-l ) * ( ONI N ) TDN( N- 1 ) ) ) *CRNI 
8= 1 1 -2. DOT AH  N) *(ON(N)-DN(N- 1 ) ) ) /( DN I N ) *DN I N— 1 ) ) +A2 I N) ) *CRN I 
1 ♦A4(N)/DS 

C=(2.0DO*A1(N)*DNIN-1  ) ) / t'ONt  N ) *(  ONI  N )4DN(  N- 1 ) > )*CRNI 
D=-(  Ml  NN(N)*A1 (N)*V1N(N) +A2(N)*K1(N)  ) * ( 1 . C DO-CRN I ) 

1 -A3(N)+A4(NI*«rl(N)/DS 

E(N)=-C/(BtA*E(N-1 ) ) 

10  F(N)x( D-A4FIN-1  1 )/( B4A4EIN-1  ) ) 

IF  (SSFAC)  11.11.12 

12  SK 1=CS 14 (ONI  IM-1  )42.0D0*DN( IM)  )/ I DN(  IM ) * I ON  I I M ) ♦ DN I I M- 1 ) ) I 

1 -(ONI  IM-1  >+DN(  1 M ) ) *EI  IM)/ IONI  IM-1  )*DN(  IM)  ) 

2 -ONI IM)*(B*E( I M ) *C  ) / ( A40NI IM-1 >*<ON( I M)+0N( IM-1 ) ) ) 

SK2=-CS2 ♦ I DN I IM-1 ) TON I IM  ) I *F I I M ) /( DNI I M- 1 ) 40NI  IM) ) 

1 TONI  I M)*(B*F( I M)-D>/( A*DN( IM-1  )* (ONI  IM-1  )+DN(  IM ) ) ) 

W2(IE)=SK2/SK1 
GO  TO  13 
11  W 21  I £ ) =END 

13  KON=IM 

DO  20  N=2* IE 

W2IK0N )=E(KON)*W2(KONt l ) +F IKON ) 

20  KON=KON-l 

W2I 1 )=E1*W2(2 ) +F1 
GO  TO  (21.5)  . I8ACX 
21  CONTINUE 

00  30  N=2.IM 

W2NNI N ) = 2.0D0*  I W2I N + 1 ) /DNI N) +W2C N- 1 ) /DNI  N- 1 ) )/( ONI N ) *0N( N- 1 ) > 

1 —2 .000  4W2(N)/( ON  I N ) *DN 1 N— 1 ) ) 

30  W2NIN)=!DN<N-1 >*W2(N+l ) /ONI N I -ONI N> * W2 ( N-l  )/0N(N-l  ) )/(DN(N)  + 

1 ONIN-1)  M-IDN(N)-DNIN-l) ) * W2I N ) /( ONI N ) *DN( N- 1 ) ) 

W2NI  1 )=-W2(  1 )* I ON  1 2 ) + 2.0DO*DN( 1 ) )/(DN( 1 > *( DN (2 ) ♦ DN ( 1 )) ) 

1 +W2(2>*(DN(2)+DN( 1) )/ I DNI 2 ) *DN( 1 ) ) 

2 -W2I3) *DN( 1 )/( DNI 2)*( ONI  1 )40N( 2) ) > 

W2NI IE  > = W2( IE)* (DNI IM-I )+2.0D0*0N( I M I ) / 1 DN I IM )* l DN I I M ) +DN(  IM-I ) ) ) 

1 —M2!  I E- 1 ) * I DNI IM-1 )♦ ON (IM) )/(DN(  I M ) * ONI  I M- 1 ) ) 

2 TW2I IE-2 ) *DN( IM)/(0N1 IM-1 )*( ONI IM) »ON( IM-1 ) ) ) 

W2NNU  ) =— W2NI  1 ) * I ON  ( 2 ) + 2.  DO* DNI  1 ) )/(ON(l  ) * I ONI  2 ) +DNI  1 ) ) ) 

1 4W2NI 2 ) * I ON  I 2 ) TON I 1)  )/(DN(2)*DN(  1) ) 

2 — W2NI 3 ) *DN( 1)/(DN(2)*I ON  1 1 1 + DNI2 ) ) ) 

W2NNI  IE) =W2N(  IE) *(DN( I M- 1) +2.  00* DN{ IM)  )/(DN(IM)*(ON(  I M ) TON ( IM-1 ) ) ) 

2 ♦ W2N I IE-2 ) *DN( I Ml/IDNI IM-1 ) * I ON I I M I +DN ( IM-1  > ) ) 

1 -W2NI IE-1 >*(DN( IM-1 )+DN I IM) )/(ON( IM)*DN( IM-1 ) ) 

GO  TO  100 
5 CONTINUE 
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0040 

0041 

0042 
00  43 

0044 

0045 


K2N(  !)=O.OODO 
00  101  N=2.IE 

W2NIN >=<W2( N» -W2 ( N-l ) )/ON<N-l ) 
*2NN<  N>={ W2N(  N) -W2N(N- 1) l/ON(N- 1 ) 
101  CONTINUE 

W2NNI 1 )=(W2N(2)-W2N( 1 ) ) /ON ( 1 ) 

100  CONTINUE 
RETURN 
END 


DATE  = 76296 


I 


05/03/31 
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0006 

0007 

0008 

00C9 

0010 
0011 
0012 
00  13 
0014 
00  15 


IV  G LEVEL  21  SOUND  DATE  = 76296  05/03/31 

SUBROUT  I NE  BOUND ( W 1 NN, W 1 K.ttl ,CB1  ,CB2 . E 1 . FI . CRN I > 

IMPLICIT  REAL  * 8 <A-H.  O-Z) 

COMMON  /PEQS/  DS  . DNt  20  1 I • IM  • I E • A 1 1 20  1 ) . A2(  201  ) . A3I  20 1 ) . A4  ( 20 1 ) . 

1 XN( 202 ) 

01  MENS  ION  W1NN<  201 1 ,W1N( 201 ) >W1 ( 201 > 

A=  (2»  0 DO -A 1(21*  DN( 2 )I/(DN<!) *(DN(2)4DN( 1 > 1 >*CRNl 

B= ( ( -2.00+ AI < 2 > * (DN{2>  -DN<1)  1 > /< ONI  2 I *DN ( 1)  I+A2I2)  » *CRN I +A 4 ( 2 ) /DS 
C~ (2 • 0 DO  + A 1 ( 2 ) 4DN( 1)1/1 DN I 2 ) *(  DN<  2 I ♦ON< 1 ) ) )*CRNI 
D=-(W1 NN(2 ) +A1  (2 >*WIN{2 I + A2  ( 2 > *W  1 < 2 I)  *< l.OOO-CRNl ) 

1 — A3< 2)4A4(2I*V1I2)/0S 

XK1  = CB1-(DN(2) 42.0DC*ON(  l ) )/( DN( 1 > * ( DN < 2 > +DN < 1 ) ) ) 4A  AON « 1 )/(C*DN( 2) 

1 «(0N(2 >+DN( 1) ) I 

XK2-I DN(2I+DN(  1 ) )/(0N( 2) *DN( 1 >)+B*DN( 1 )/ < C*DN( 2 ) * ( ON ( 2 1 +DN ( 1 ))  > 
XK3=CB2-D*ON(l >/ (C*DN< 2 ) ♦ < DN ( 2 >+DN C 1 ) ) ) 

El  =-  XK2/XK 1 
Fl=t-XK3/XK1 
RETURN 
END 
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0030 


IV  G LEVEL  21  GEOM  DATE  * 76296 

SUBROUTINE  GEOM < S . OS, RS . CK ,CSF . S l F , XB ) 

C SPHERE-CONE 

IMPLICIT  REAL*  8 < A-H,  O- Z I 

COMMON/BASU/  XS3  , CONE 

COMMON/  PUSHY/1  DERIV1 .THMAX 

F*12. 5000 
XS=S+DS 
XS3=XS 

SMAX=3.  141 5926535  897932  DO /2.  0D0  — THMAX 
I FIXS.GT. SMAX ) GO  TO  100 

RS=DSI N( XS1 

CK=1 • 0D0— 1 • OOO/I  1 • 000 + 0 E XP ( — FA ( XS-SMAXI  ) ) 

CK=l.DO 
CCNE— — 1 « ODO 
CSF=rs 

S*.F=OCOS(XS) 

XB  = 1. OOO-OCOSI  XSJ 
GO  TO  200 
100  CONTINUE 

X S 1 = X S—  SM  A X 
CSF=DCOS(THMAX ) 

SIF=OS  INI THMAX) 

CK= 1.00 0-1. 000/(1. ODOFOEXP{-F*{ XS-SMAX) ) ) 

CK=0.  ODO 
CONE= 1. OOO 

RS=OS I N(SMAX) +XSI *DSIN( THMAX) 

XB=1. 00  0— OCOS ( SMAX) FXS1 *DCOSI THMAX) 

300  FORMAT <4F10.  5 > 

200  CONTINUE 

RETURN 
END 


05/03/31 
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0005 

0006 

0007 

0008 

0009 

0010 
00  1 1 
00  12 
0013 
00  14 

00  15 
0016 

0017 

0018 

0019 

0020 
0021 
0022 
0023 
00  24 

0025 

0026 
0027 
0928 

0029 

0030 

0031 
00  32 

0033 

0034 

0035 

0036 
0 0 37 

0038 

0039 

0040 

0041 

0042 

0043 

0044 


IV  G LEVEL  2 I 


SHVALS 


DATE 


76296 


05/03/31 


SUBROUTINE  SHVALS  ( SP • CP.  SP8.  CPB.  TTSH,  VRSH.  URSH,  PPSH,  ID* 
IMPLICIT  REAL48  <A-H.  O-ZI 


COMMON  / INSH/  CONO 
1 EPS 

COMMON/OUTSH/  PPS 


PPSl 

PSP 

PPS2 


GAM 

RMAC 

RPS 

RRSl 

RRS2 

RSP 


UPSH 

VI  SCO 

UUS1 

UUS2 

USP 

uus 


XNS 


VVS 

VVS1 

VVS2 

VSP 


s 

TP  SH 
TTS 
TTSI 
TSP 
TTS2 

DS > ONI  20 1 ), IM. IE. All  201) , A2( 201 ) .A3(  201  * .A4I201  I . 

XN1202* 

= GAM  4 1 • 0D0 
= GAM  - 1 .000 
* RMAC  * RMAC 
- 4. ODO/C  GAMP* GAMP  * 

= SP  * SP 
» EPS  * EPS 
= RMACO  * RMACO  4 SPQ 
SP  * CP  / (SP  ♦ EPSO  * VISCO  * UPSH/XNS) 

< (URSH-CP)**2  4 F CGQ* GA  M4  SPQ  4(2*  ODO/GAMM— F0GQ4GAMM J/RMACQ 
- FOGQ/DEN*  *0.5D0*SP/(SP4EPSQ*C0N0*TPSH/XNS I 
1 2 .0004  SPQ  - GAMM/ (GAM4RMAC0) ► / GAMP 


10 


COMMON  /PEQS/ 

t 

GAMP 
GAMM 
RMACO 
FOGQ 
SPQ 
EPSQ 
OEN 
URSH  a 
TTSH  > 

1 

PPSH  a 
RR  SH  * GAM  * PPSH 
VRSH  = -SP  / RRSH 
GO  TO  (20.51  . ID 
CONTINUE 
TTS2  «*  TTSH 
PPS2  = PPSH 
RRS2  = RRSH 
UUS2  = URSH  * 
VVS2  aa  -URSH  * 

IF  (S  . GE.  .0001  ) 
UUS1  aa  -UUS 2 
VVS1  = VVS2 
TTSI  = TTS2 
PPSl  as  PPS  2 
RRSl  sa  RRS2 
CONTINUE 


/ (GAMM  4 TTSH* 


SPB  4 VRSH  • 
CPS  4 VRSH  4 
GO  TO  10 


CPB 

SPB 


UUS 
VVS 
TTS 
PPS 
RRS 
USP 
VSP 
TSP 
PSP 
RSP 
20  CONTINUE 
RETURN 
END 


(UUS1 
(VVS1 
(TTSI 
(PPSl 
(RRSl 
( UUS2 
( V VS2 
(TTS2 
( PPS2 
(RRS2 


UUS2  t 
VVS2  » 
TTS2  » 
PPS2 ) 
RRS2  > 
UUS1  ) 
VVS1  I 
TTSI  ) 
PPSl  * 
RRSl  ) 


2a 

2. 

2. 

2. 

2a 

DS 

DS 

DS 

DS 

DS 


0D0 

000 

0D0 

0D0 

ODO 
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0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 
0000 
00  10 
00  1 1 


IV  G LEVEL  21  BOUND1  OATE  = 76296  05/03/31 

SUBROUTINE  30UN0 1 ( EE  1 .FF1 ) 

IMPLICIT  REAL *8  (A-H,  O-Z) 

COMMON  /PE QS/  OS . ONI  20 1 ) . I M . I E . A 1 ( 20 1 ) , A2 ( 20 1 I . A3 < 20 1 > . A4 ( 20 1 ) . 

1 XN1202) 

A =1.0  DO 

B=-2.  000*DS*DS*A2C2I 
C=1.0D0 

D=-DS*DS*A3<  2 ) 

EE1  = < B+4.0D04C l/( 3. 000* C- A) 

FF1=-D/(3.0D0*C-A> 

RETURN 

END 
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0001 

0002 

0003 

0004 

0005 

0006 
0007 
0003 

0009 

0010 
0011 
0012 
00  13 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 
0029 
00  30 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 


SUBROUTINE  PUSHPA (YNSH,DS  ,IENP,YNSP> 

IMPLICIT  REAL*8  <A-H,  0-Z) 

COMMCN/KINNI/  XNSHI 110) f XNSP I 1 10 ) . XNSPP II 10 ) 

COMMON/PUSHY/  OER I VI ,THMAX 
COMMON/CON/  NJNC,NJ1,RUMP 
DIMENSION  YNSH(llO) »YNSP(110) 

DIMENSION  ACKP(llO)  .APHIUOJ 

COMMON/BASU/  XS3  , CONE 

CCMMON/MANIS/  AXSH 1 1 10 ) , AXS P 1 1 10 ) , AXSPP 1 1 10 ) 

XS=O.ODO 

I END1 = I EN0*1 

DO  500  I =2  » I END1 

CALL  GEOM  I XS ,DS, RS , CK ,CSF , SI F , XB ) 

XNSHI I ) = l YNSH l I J-RSI/CSF 
A XSH { 1 ) =XB— XNSH ( I I *SI  F 
ACKPII 1 = 11 .COO+CK*XNSHI I ) ) 

APHII )=DAPCOS(CSF) 

XS=XS3 

500  CONTINUE 

XNSHI 1) =<4.000*XNSH12)-XNSH< 3M/3.0D0 
AXSHI 1) =- XNSHI 1 ) 

DC  12  N=2»  I END 

IF  I PUMP. LT. 0.000)  GO  TO  510 

IF(N.EQ.NJl)  GO  TO  501 
IFtN.EQ.NJNC)  GO  TO  504 

510  CONTINUE 

AXSP(N)=( AXSHI Ntl J-A XSH IN-1) )/(2.0D0*0S» 

GO  TO  502 

501  CONTINUE 

AXSP[NJ1)=(3.0DO*AXSH(NJ1)-4.0DO*AXSHINJ1-1) +AXSH ( NJ1 -2 >> / 
1I2.0D0*DS) 

GO  TO  503 

504  CONTINUE 

AXSPINJNC) =(4.0DO*AXSH(NJNC+1)-AXSH{NJMC*2)-3.0DO*AXSHINJNC))/ 

II 2.000*DS) 

503  CONTINUE 

502  CONTINUE 
12  CONTINUE 

XNS  P 1 1 ) =0 .ODO 
AXSP I 1 ) =0. OD  0 

AXSPI IEND1)=I3.0D0*AXSH(IEND1)-4.0D0*AXSH{ IEN01-1 )+AXSH(I END1-2I) 
l/I2.0D0*DS) 

DO  860  I = 1 » I END1 

IFII.EQ.l)  GO  TO  81 
TALP=YMSP(II/AXSPII) 

XNSPII )*ACKP(I 1*1  I TALP-DT  AN  I APH I I 1 ) ) / 1 1 .O+TAL P*DT AN  I APH 11 ) ) I) 

GO  TO  511 

512  CONTINUE 

IF! I.LT-IE  NO  1 1 XNSP (I  I = I XNSH I I + 1 ) -XNSHI t-l ) ) / I DS*2 .ODO ) 

IF! I .EQ. IENDI ) XNSP I I)  = I 3. ODO* XNSHI I END  1 ) -4. ODO* XNSH I I END 1-1 ) 
l+XNSHI IEND1-2) )/<2.0D0*DSJ 

511  CONTINUE 
81  CONTINUE 
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860  CONTINUE 

XNSPP ( 1)  = (2.000*XNSHI  1 ) -5  .OODO*XNS  H ( 2 ) +9  .ODO*XNSH  ( 31  -XNSH  I 9 1 I 
1/<DS*C$» 

RETURN 

ENO 
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DERI  V 


DATE  = 76296  05/03/31 


0001 

0002 

0003 

0904 

0005 

0006 

0007 

0008 

0009 

0010 
00  1 1 
00  12 
00  1 3 
00  14 
0015 
00  16 
00  17 
0018 
00  19 
0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 
0029 
00  30 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 
00  44 

0045 

0046 

0047 

0048 
0040 

00  50 

0051 


SUBROUTINE  DERIVI DS  » I END • I E ND  1 . AXNSH  , A XNSP  . A XNSPP  ) 

IMPLICIT  REAL  *8  (A-H,  0- Z ) 

D IMENSI ON  XI SP( 1 1 0) ,A 1SP( 1 10) , XI  SHI  1 10) 

DIMENSION  YNSH<  HOI  .YNSP<  1 1 0 1 , YNSPP  { 1 1 0 ) 

DIMENSION  ACKPdlOl  .APH(llO) 

COMMON/CON/  NJNC.NJ1 .RUMP 
COMMON/BASU/  XS3  . CONE 

COMMON/KI NN 1 / X NSH ( 1 1 0 > . XNSP (110). XNSP  P(UO) 

COMMON/MANIS/  AXSHt 110) »AXSP( 1101 . AXSPP (110) 
COMMON/PUSHY/  DER1V1 .THMAX 

DIMENSION  AXNSH t 1 19) . AXNSP( 119) . A XNSPP ( 119) 

DER IV 1— 1 .ODD 
AHALF=1  .ODO 
A HALF  =— 1 • ODO 
IE  ND1  = I ENDA1 


READ( 5.80) 

< XNSH( I ) 

1=1.8) 

READ! 5. 80 ) 

< XNSH ( I ) 

1=9.16) 

READ! 5,81 ) 

< XNSH( I ) 

1=17,221 

READ! 5.80) 

( XNSH( I ) 

1=23.30) 

READ(S.aO) 

< XNSH ( I ) 

1=31. 38) 

READ ( 5.81 ) 

(XNSH ( I ) 

1=39.44) 

READ( 5.80) 

( XNSH ( I ) 

1=45.52) 

READ ( 5.80 ) 

( XNSH ( I ) 

1=53. 60) 

READ ( 5.82) 

XNSH(61 ) 

82  FORM  AT ( F 1 0 • 6 ) 

80  F0RMAT(8F1 0. 6 ) 

81  FORMAT ( 6F 10.6) 

XS=0. ODO 

DO  500  I =2  * I END! 

CALL  GEOM  C XS , DS , RS . CK ,C SF , S IF , XB ) 

YNSH< I )=RS+XNSH( I >«CSF 
A XSH(  I )=XB-XNSH(I  >*SIF 
ACKP( I )=(1.0DO+CK*XNSH{ I ) ) 

A FH  ( I ) = DARCOS(CSF) 

I F( I • EO • 2 ) YS1=RS 

XS=XS3 

500  CONTINUE 

YNSH< 1)=0.0D0 
AXSHC  1)=-XNSH(1  ) 

DO  12  N=2.IEN0 

lF(RUMP.LT.O.ODO)  GO  TO  510 

IFIN.EO.NJll  GO  TO  501 

IF(N. EQ.NJNC)  GO  TO  504 

510  CONTINUE 

YNSPI N)  =( YNSHlNtl ) -YNSH(N-1 ) )/<  2.0D0*0S ) 

AXSP(N)=(AXSH(N+1 )-AXSH(N-l ) )/<2.0D0*DS) 

GO  TO  502 
501  CONTINUE 

AXSP<  NJ1 )=(3.0D0*AXSH( NJ1 )-4. OD 0* AXSH( N Jl- 1 )♦ AXSH ( NJ 1 -2 ) ) / 
1 (2.0D04DS) 

YNSP( NJ1 )=<3.0DO*YNSH{ NJ1 )- 4.  00  0*YNSH(  N J l- 1 ) ♦ YNSH ( N J 1-2 > ) / 

1 (2 .ODO  *DS) 

GO  TO  50  3 
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0052 

0053 

OC  54 

0055 

0056 

0057 
0 0 58 
0059 
00  60 
00  6 1 
00  62 
00  63 

00  64 

0065 

0066 
00  67 

00  68 

0069 

0070 

0071 

0072 

0073 

0074 

0075 
0 0 76 

0077 
00  78 
0 0 79 
00  80 
0081 
0082 

0083 

0084 
0 0 85 
0086 

0087 

0088 

0089 

0090 

0091 

0092 


IV  G LEVEL  21  DERIV  OATE  = 76296  05/03/31 

504  CONTINUE 

AXSPI NJNCI=(4.0DO * AX SHI NJNC+1 I - AXSH I NJNC  + 2 > -3. OOO * A XSH ( N J NC ) )/ 

1 <2 .0D0*D5) 

YNSPI NJNC  )=l 4. ODO* YNSH I NJNC  + 1 ) - YNSH I N JNC+2 I -3 . ODO * YNSH I NJ NC  > > / 

1 12.0D0*DS> 

503  CONTINUE 

502  CONTINUE 
12  CONTINUE 

YNSPI 1 1 = 14. 0D0*  YNSH ( 2 ) -YNSH! 3)— 3.0D0*Y  NS HI  1 ) ) / ( 2 . ODO  *DS ) 

YNSPI  1 ) — XNSHI  1 l+YSl/OS 
YNSPI 1 >=XNSHI  1 )+l «ODO 
XNSPI  1) =0.000 
AXSPI 11=0.000 

AXSPI  IEND1  )=I3.0D0*AXSHUEND1  > -4 .0 DO * AX SH I I END  1-1  l + AXSHI  I END  1-2)1 
l/l 2.0D04DS) 

YNSPI  I E NO  1 1 = 1 3.0D0+YNSHI IEND1 ) - 4. 0 DO * YNSH I I END  1-1  l + YNSHI I END 1-2) ) 
1/t  2.0D0*DS) 

DO  700  1 = 1.  1 END  1 

IFII.EQ.l)  GO  TO  600 

TALP= YNSPI  I > /AXSPI I ) 

C IF1 RUMP.LT .0.000 ) GO  TO  512 

XNSPI I ) =ACKP!  I ) *1  I TALP-DTANI APH I I ) ) I / I 1 . 0 + T ALP *DT AN I APH I 1)1)) 

GO  TO  511 

512  CONTINUE 

IFI  I. LT. I END  1 ) XNSPI 1 1 =! XNSHI 1+1  )- XNSHI I ) )/ I 2. 0D0*DS ) 

IF(  I .EO. I END  1 ) XNSPI I )= I 3. ODO* XNSHI  I END  1 )-4. ODO* XNSHI  IEND1-1 ) 
1+XNSHI  IEND1-2I  )/I2.0DO*DS) 

511  CONTINUE 

600  CONTINUE 
700  CONTINUE 

DO  28  N=2  * IE NO 

C IFI RUMP.LT. 0.000  I GO  TO  530 

IFIN.EO.NJ1)  GO  TO  S3 1 

IFIN.EQ.NJNC)  GO  TO  532 

530  CONTINUE 

YNSPP IN)  = IYNSHIN+1  ) + Y NS  HI N— 1 ) -2 . 0 DO * YNS H I N » >/IDS*OS) 

GO  TO  533 

531  CONTINUE 

YNSPP IN > = < 2. ODO* YNSH! N) +4. ODO* YNSH! N-2) -5. ODO* YNSH IN- 1 ) -YNSH!  N- 3) 

1 )/  I DS*  OS  ) 

GO  TO  533 

532  CONTINUE 

YNSPPlN)=I2.3D0*YNSHtN)-5.0D0*YNSH{N+l ) +4 . ODO* YNSH I N+ 2 ) -Y NSHI N+ 3 ) 

1 )/  IDS*  DS  ) 

533  CONTINUE 

28  CONTINUE 

YNSPP 11 ) = I 4. 0 DO* YNSPI 2 ) -YNSPI 3) -3. ODO *Y NS PI  1))/I2.0DO*DS) 

YNSPPI1 )=0.000 

YNSPP! I END1  ) = I 3 • 0 DO  *YNS  P I I END!  > -4 .0  DO *Y NSPI  I END  1— 1 )♦ YNSPI  I END  1-2) 

1 ) / 1 2. 0 D0*DS ) 

YNSPP!  IEN01  ) = ! 2.0DO*YNSHI I END  1 )- 5. ODO* YN SH I IEND1-1  ) + 4. OOO * YNSH I 
HENOl -2  I -YNSH!  I END1  — 3 ) ) / I DS*DS  ) 

C *+*♦*+**»♦**:***«**  + **  *********************  *********** 
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0093 

0094 

0095 

0096 

0097 

00  96 

0099 

0100 
Old 
0102 

0103 

0104 

0105 

0106 
0107 

01  OB 
01  09 
0110 
0111 
01  12 
0113 
0 114 

0115 

0116 
0 117 
0118 
0119 
0 120 
0 12  1 
0 1 22 
0123 
Cl  24 

0125 

0126 
0127 
0 1 28 
0129 
0 130 
0131 
0 l 32 

0 1 33 

0134 

0135 

0136 

01  37 

0138 

0139 

0140 

0141 

0 142 

0143 

0144 


I F< AHALF.LT.O.ODO I GO  TO  6006 

APARA=-1»ODO 

APARA=1 .ODO 

6019  CONTINUE 

DO  6009  1 = 1.1  END 

AXNSHU+I  ) = [YNSH(  I l+YNSHd+ll  1/2.  ODO 
A XNSP ( 1*1 1=(YNSP(  I l+YNSPI 1 + 1 1 1/2.000 
A XNSPP ( 1+ I 1 = ( YNSPP ( I M-YNSPP  ( 1 + 1 I 1/2.0  00 
A 1SP(  l+ll=( AXSP( I 1 + AXSP  <1  + 1 1 1/2. ODO 
X1SHI I* I !=<XNSH< I l+XNSHI 1+1 1 1/2. ODO 
X1SPC  I+I1=IXNSP< I I +XNSP ( 1 + 1 I 1/2. ODO 
6009  CONTINUE 

NJJi=NJ 1+NJ1 
I 2END  = 2*IEND«-  1 
I I = i 

DO  6007  1=1 . I2END .2 

A XNSH< I l = YNSH (II) 

A XNSP  < I 1 = YNSP (III 
A XNSPP ( II=YNSPP( I I 1 
X 1SH(  I } = XNSH<  I I 1 
X 1SP(  I )=XNSP< 1 I 1 
A 1 SP(  I 1=AXSP( I I I 
I 1=1  I +1 

6007  CONTINUE 
N J1  = N JJ  1 
NJNC=NJJ1*1 
I END1 =1 2END 
I END= I END  1—1 
DO  6011  N= 1*1  END  1 

XNSH ( N 1 = X1 SH ( N I 
A XSP  ( N 1=A  1 SP(  N I 
XNSP<  N!=X1SP<  N 1 
6011  CONTINUE 

AXNSHC NJ J1 ) = ( 2. ODO* AXNSH(NJJ 1-1 1 — AXNSH (NJJ1-2) ) 

A XNSP ( N J J 1 ! = < 2. ODO* AXNSP(NJJ1-1  1- AXNSP ( NJ J 1 -2 1 1 
A XNSPP ( NJJ1 1 = ( 2.  000*AXN5PP( NJJl-t 1- AXNSPP ( N J J 1-2  I I 
X NSH(  NJ J1 1=<  2.0D0*XNSH(  NJJ1-1  l-XNSHINJJ 1-211 
XNSP{ NJ  J 1 1=(2. ODO*XNSP( NJJ1-1  1 -X NSP ( NJJ 1 -2 1 1 
AXSP(NJJ11=(2.0D0*AXSP( NJJ1-1 1- AXSP ( N JJ 1-2 1 1 
IF< APARA.GT.0.0DC1  GO  TO  6020 
DO  6021  1=1.1  END  l 

YNSH ( I l=AXNSH(  I I 
YNSPt  I 1=AXNSP( I 1 
YNSPP ( I l = AXNSPP ( I 1 
6021  CONTINUE 

GO  TO  6019 

6020  CONTINUE. 

6006  CONTINUE 

IF< AHALF.GT.O .ODO 1 GO  TO  6010 

C ***************************************************** 

DO  1200  N= 1 » I END  1 

RE AO (5.1 201  I Y NSH ( N 1 ,YNSP(N1 .YNSPP ( N1  , AXSPf  N! , XNSH( N 1 ,XNSP( N1 

1200  CONTINUE 


169 


AEDC-TR-77-20 


FORTRAN  IV 

0145 

0146 

0147 

0148 

0149 
01  50 
0151 
Cl  52 
0 1 53 
0 1 54 
0 1 55 


G LEVEL  21  DERIV  DATE  = 76296  05/03/31 

XNSPPI 1 )={ 2» 0D0*XNSH( 1 ) -5. 0 000* XNSHI 2 ) +4 . 0D0*XNSH ( 3 )-XNSH ( 4 ) I 
l/( DS*DS) 

00  6008  I =1.1 END 1 

AXNSHI I )=YNSH( l ) 

A XNSP { I l = YNSP< 1 ) 

A XNSPPI I l±YNSPP( l » 

6008  CONTINUE 

6010  CONTINUE 
1201  FORMATI6F12.8) 

ANSH= YNSHt IEN01 ) 

RETURN 

END 
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C0C1 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
00  1 1 
00  12 
00  13 
00  14 

0015 

0016 
0017 
00  18 

0019 

0020 
0021 
0022 
0023 
00  24 

0025 

0026 

0027 

0028 
0029 
00  30 

0031 

0032 

0033 

0034 

0035 

0036 
0 0 37 

0038 

0039 
00  40 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 
0 0 50 


IV  G LEVEL  21  MANISH  DATE  = 76296  05/03/31 

SUBROUT INE  MANtSHt ARXX  , ARX , AR , RL . DX , AAK3 . A AKA , A AK5 , AM) 

IMPLICIT  REAL*8  ( A-H,  0-2) 

COMMON/MANU/EE 1 ,FF1 . IENO, I END  1 . AAAI  C 1 1 0 ) , AA A2 < 11 0 ) . A A A3<  110) 

1 .AAA4< 110) 

COMMON/CON/  NJNC .NJ1 , RUMP 
DIMENSION  ARXX ( 110). ARX(110),AR <1101 

DIMENSION  E ( 120  I ,F{ 1 20 ) 

EUI-EE1 
F( 1 )=FF1 
XI=0.  50 DO 
R0=0.  0D0 
I END1  — I END  4-1 
IM=IEND1-1 

IF(RUMP.GT.O.ODO)  GO  TO  501 
AM=1 . ODO 
AAK3=1. ODO 
AAK1=1.000 
AAK4=0.0D0 
AAK5=0. 000 
501  CONTINUE 

DO  10  N=2  * IM 

ALP  1 B = A A A 1 { N ) 

ALP1 A=A AA 1 ( N ) 

ALP2= AAA2 ( N) 

ALP3B=AAA3(N) 

ALP3 A — A AA3  < N ) 

AK1=1 .000 
AK3=1 .000 
A K4=0 . 0 DO 
A K5=0 • ODO 

IF(N.LT.NJl)  GO  TO  11 

I F(N. GT.NJNC)  GO  TO  11 

AK1=AM 

AK3=  AAK3 

AK4=A AKA 

AK5=AAK5 

IF(N.EQ.NJl)  .GO'  TO  12 

AP1=XI*( AK1-DX4XI *AK5/2 .ODO )♦(  1 .ODO-X I ) 

AP2=X I* ( 1 . ODO-X I ) *< AX1-CXAXI* AK5/2. ODO) ♦ X I * X I * AK3/2. 0 DO ♦ 

1 ( 1 .ODO-X I ) **2.  0 00/2.000 
AP3=A  PI +2, 000  * AP2 

A=2.0D0/( DX*DX*AP3)-ALPIB/(DX*AP3 I 

B=-2.  COO* ( 1 . ODO + API  >/t  0X*DX*AP3  1+ALP1 B* ( 1 .0D0-2. ODO * AP2 > / ( OX* AP 3 ) 

1 4ALP2 

C=APl*2.0D0/(  DX*DX*AP3) *2.000 * ALP 1 B*AP2/ ( DX*AP3 ) 

D=- ALP3 B+  X I *X I * AK4/AP3  — OX* ALP1 B*X I *X I *AK4/ ( 2.0DO*AP3) 

GO  TO  13 
12  CONTINUE 

A S=( 1 . ODO  — XI )*(1.0D0-DX*( 1.000- XI ) * AK5/ < 2 . 0 DO* AK3 ) )/AKl 
AP1=XI*AS 

AP2=X I*XI/2. ODO*XI*AS*( 1. ODO-X I )**2.0D0/<  2.0D0+AK3) 

AP3=AP1*2.0DO* AP2 

A=APl*2,0D0/{  DX*DX*AP3)-2.0D0*ALP1A*AP2/<DX*AP3) 
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fortran 

0051 

0 0 52 
0053 
0 0 5* 

0055 

0056 
0 0 57 
0058 
0 0 59 
0 0 60 
0061 
0062 

0063 

0064 

0065 
0 0 66 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 
0081 
0082 
0 0 83 

0084 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 


IV  G LEVEL  21  MANISH  DATE  = 76296  05/03/31 

B=-2.000*< 1.0D0+AP1 » / ( DX*0X*AP3 ) - ALP1 A* < 1 . 000-2. ODO* AP2 I / < DX* AP3 1 
l*ALP2 

C=2.  ODO/(DX*DX*AP3)*ALP1A/{DX*AP31 
A14=AK4/AK3 

D=-ALP3 A-ALP1  A*OX*(C1.0DO-XI)**2.0DO)*A14/(2.0DO*AP3I-(1.0DO-XI) 
1**2.0D0*A14/AP3 
GO  TO  14 
11  CONTINUE 

A = 1.0DO/<  DXOXl-AAAH  N)  /{2.  0D0*DX> 

B=-2. ODO/ ( DX  * DX I +AAA2IN I 
C=1.0D0/IDX*DX)*AAA1 (N)/(2.  0D0*DX) 

0=-AAA31N» 

14  CONTINUE 

13  CONTINUE 

E <N)=-C/( B4A*ECN-1  ) 1 
F(N>=CD-A*F(N-1 >)/<e+A*E(N-l ) ) 

10  CONTINUE 
KON= I M 
AR< I END1  i =RL 
AR< 1 »=R0 

DO  20  N=2.IEND1 
A RC  KONI =E { KON  1 * AR  ( KON  4 I )4F(K0N> 

20  KON=  KON—  l 

C CALCULATION  OF  DERIVATIVES 
DO  30  N = 2 . I M 
IF(N.EQ.NJl)  GO  TO  31 

IF(N. EQ.NJNCI  GO  TO  40 

A RX ( N 1 = (AR(N* 1 ) — AR(N-1 ) l/I2.0D0*DX) 

ARXXI N)=<  ARIN41  )*AR{N-  1 1-2. 00 0* AR ( N ) )/( DX*DX) 

GO  TO  50 
31  CONTINUE 

AK1=AM 
A K3  = A AK3 
AK4-AAK4 
AX5=AAK5 

AS=t  1 .ODO-XI 1 *( 1 .0D0-DX*< 1.0D0-XI 1* AK5/ ( 2 . ODO* AK3 ) ) /AK 1 
API =X 1+ AS 

AP2=X I*XI/2.0D0*XI*AS*(  1. ODO-XI )**2.0DO/<2. 0D0*AK3) 

AP3=AP1+2.0D0*AP2 

A14=AK4/AK3 

ARX(N)=(AR(N*1 1-2.0C0* ARCN-I 1 *AP2-AR ( N 1 *( 1 . 000-2. ODO* AP2 1 +DX*DX 
1*(  1. ODO-XI >*C 1 . ODO-XI 1 *A 14  >/(DX* AP3I 
ARXXlN  > = < AP1*AR(N— 1 I + AR ( N* 1 1 -AR ( N > * ( AP 1 ♦ 1 . 0 00 1 +DX*DX* ( l. ODO-XI  1* 
1I1.0D0  — XIl*A14/2.0D0)/( DX*DX  *AP3/2 .0  DO  1 
GO  TO  SO 
40  CONTINUE 
A K1 *A  M 
AK3-AAK3 
AK4=AAK4 
AK5=AAKS 

AP1=X  I* (AK1-DX*XI * AK5/2.0D0 l + < 1 .ODO-XI I 

AP2  = X I* (1 . ODO-XI  1 * { AK 1 — CX  *X I * AK  5/2 • ODO )+  XI*XI*AK3/2. 0D0  + 

1< 1 .ODO-XI 1**2. 0D0/2.CD0 
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FORTRAN  IV 

OOS0 

0099 

0100 

0101 

0102 

0103 

0109 

0105 

0106 

0107 

0108 
0109 


G LEVEL  21  MANISH  DATE  = 76299  23/19/20 

AP3=AP1+2.0C0*AP2 

ARX(N)=(2.0D0*AP2*AR(N+1)+AR(N>*(  1.000-2. 000*AP2 1 -AR ( N-l) »OX*OX 
1*XI»XI*AK9/2.C'C0)/(0X«AP3) 

ARXX(N)  = (AR(N-1)+AR(N+1)*AP1-AR(N)*(1.0D0*AP1  )-OX*DX*X!*X  I *4X9 
1/2.000) /(DX*DX*AP 3/2. 000) 

50  CONTINUE 

30  CONTINUE 

AP.Xt  1>=(-3.0P0*ARI1)  + 9.0D0*AR(2)-AR13)  J/12.0D0*DX> 
ARX<IENDI)=(3.0D0*AR(IEND1)-9.0D0*AR( I END  1- 1 ) +AR ( ! ENO 1-2 ) )/ 

1 ( 2.000*DX) 

AR  XX l 1 ) =0.000 

ARXX(l)  = <2.ODO*AR(l)-5.OP0*AR(2)+9.OD0*ARm-ARI9)  )/(DX*DX> 

ARXX( IEND1)=C2.0D0*AR( I END  1 ) -5 . OCO* AR II END1-1 ) +9 .000* AR< I END1-2 ) 
1-AR( IEN01-3I )/(DX*DX) 

RETURN 

END 
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S 
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X 
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HEAT 
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CCF 
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0.075557 


CDTOT 

1.756005 

XSH 

-0.055672 


USH 
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USP 

0.772023 


VSH 

-0.135769 

VSP 

0.169557 


TSH 

0.555626 

TSP 

-0.239527 


RSH 

5.817957 

RSP 

— 0. CS7650 


PSH 

0. 750780 


PSP 

-0.510606 


CF 

0.081502 
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HEAT 
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STAN 

0. 158787 


0. 156709 


0. 521535 


CDF 

0.019760 

NSH 

0.153358 


CCP 

1.653118 

NSHP 

0.106525 


CDTOT 

1.672878 

XSH 

0.015859 


USH 

0.538613 

USP 

0.691323 


VSH 

-0.107593 

VSP 

0.236558 


TSH 

0.510673 

TSP 

-0.278653 


P.SH 

5.802299 

RSP 

-0.133720 


PSH 

0. 680865 
PSP 

-0.577657 


RSH 

0.0 


NO  ITER 
22 


PWALL 

0.921095 

RSH 

.9.155556 


PW/PO 

1.000000 

MO  ITER 
22 


PWALL 

0.899922 


PW/PO 

0.977015 


RSH  NO  ITER 
0.306586  28 


PWALL 

0.855565 

RSH 

0.556222 


PW/PO 

0.918000 


HO  ITER 
31 


PWALL 

0.761815 

RSH 
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PW/PO 

0.827075 

NO  ITER 
33 
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CF 

HEAT 

STAN 

CDF 

CDP 

COTOT 

PW  ALL 

PW/PO 

0.097430 

0.066936 

0. 140674 

0.031548 

1.547861 

1.579403 

0.654892 

0.710993 

s 

X 

0 

NSH 

NSHP 

XSH 

P.SH 

NO  ITER 

0.685665 

0.226002 

0.633188 

0.171189 

0. 153678 

0.093502 

0.741583 

33 

USH 

VS  H 

TSH 

RSH 

PSH 

0.526864 

-0.070203  0.371583 

5.781352 

0.613854 

USP 

VSP 

TSP 

RSP 

PSP 

0.607613 

0.303819  -0.299625 

-0.175594 

0.513623 

CF 

HEAT 

STAN 

CDF 

COP  , 

COTOT 

PWAU 

PW/PO 

104987  0 

.057308 

0. 120440 

0.044853  1.429514 

1.474367 

0.538326 

0.584442 

s 

X 

R 

NSH 

NS  HP 

XSH 

RSH 

NO  ITER 

622798  0 

.319827 

0.733052 

0.196135  0.210140 

0.186422 

0.876829 

32 

USH 

VSH 

TSH 

RSH 

PSH 

0.601265 

-0.022700  0.331978 

5.755171 

0.545963 

USP 

VSP 

TSP 

RSP 

PSP 

0.490385 

0.383401  -0.236747 

-0.210476 

0.491541 

CF 

HEAT 

STAN 

CDF 

CDP 

COTOT 

PW  ALL 

PW/PO 

0. 104886 

0.047424 

0.099668 

0.058590 

1.309062 

1.367651 

3.423802 

0.460106 

S 

X 

R 

NSH 

NSHP 

XSH 

RSH 

NO  ITER 

0.959931 

0.426424 

0.819152 

0.230035 

0.285009 

0.294452 

1.007627 

31 

USH 

VSH 

TSH 

RSH 

PSH 

0.660195 

0.035136  0.294536  5.724023  0 

.481780 

USP 

VSP 

TSP 

RSP 

PSP 

0.371903 

0.459203  -0.261860  -0.244099  -0 

.448872 

CF 

HEAT 

STAN 

CDF 

CDP 

CDTOT 

PW  ALL 

PW/PO 

0.098341 

0.038033 

0.079931 

0.071925 

1.196499 

1.268424 

0.32177  2 

0.349336 

s 

X 

R 

NSH 

NSHP 

XSH 

RSH 

NO  ITER 

1.09  7064 

0. 5437B9 

0.889872 

0.276469 

0.391476 

0.417661 

1.135894 

60 

USH 

VSH 

TSH 

RSH 

PSH 

0.704822 

0. 101709 

0.259543 

5.686685 

0.421797 

USP 

VSP 

TSP 

RSP 

PSP 

0.263302 

0.520587 

-0.238522 

-0.286226 

-0.408860 

CF 

HEAT 

STAN 

CDF 

COP 

CDTOT 

PWAU 

PW/PO 

0.088185  0 

1.029434 

0. 061860 

0.084468 

1.099116 

1. 183583 

0.233792 

0.253320 

S 

X 

R 

NSH 

NSHP 

XSH 

RSH 

MO  ITER 

1.234197  C 

1.669721 

0.943883 

0.343699 

0.509031 

0.556204 

1.268295 

55 

USH 

VSH 

TSH 

RSH 

PSH 

0.738785 

0.173704  0.225124  5.638285  0 

.362799 

USP 

VSP 

TSP 

RSP 

PSP 

0.226261 
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CF 

HEAT 
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CDP 

COTOT 

PWAU 

PW/PO 
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S 

X 

R 
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RSH 
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46 
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USH 

VSH 

TSH 

RSH 

PSH 

0*761950 

0.249525 

0.191662 

5.574810 

0.305444 

USP 

VSP 

TSP 

RSP 

PSP 

0.187022 
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CF 
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0.062026 

0.015841 
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S 

X 

R 
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RSH 

NO  ITER 
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36 

USH 

VSH 

TSH 

PSH 

PSH 
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CF 
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PWALL 

PW/PO 
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1.055018 

0.  080785 

0.087705 

S 

X 

R 

NSH 

NSHP 

XSH 

RSH 

MO  ITER 

1.645596 

1.073414 

1.018294 

0.596364 

0.470475 

0.995573 

1.609556 

34 

USH 

VSH 

TSH 

RSH 

PSH 

0.807850 

0.265884 

0.152159 

5.467445 

0.237738 

USP 

VSP 

TSP 

RSP 

PSP 

0.158379 

-0.106792 

-0.099747 

-0.347341 

-0.170948 

CF 

0.040721 

HEAT 

0.010799 

STAN 

0.022696 

COF 

0. 125728 

COP 

0.909735 

COTOT 

1.035462 

PWALL 

0.071915 

PW/PO 

0.078075 

S 

1.762729 

X 

1 . 209373 

R 

1.036193 

NSH 

0.658261 

NSHP 

0.432264 

XSH 

1. 123453 

RSH 

1.688823 

NO  ITER 
43 

USH 

VSH 

TSH 

RSH 

PSH 

0.827629 

0.251828 

0.139646 

5.420094 

0.216294 

USP 

VSP 

TSP 

RSP 

PSP 

0. 127609 

-0.097378 

-0.081236 

-0.335642 

-0.139215 

CF 

U. 036198 

HE  4T 

0.009647 

STAN 

0.020275 

COF 

0.131422 

COP 

0.883284 

CDTOT 

1.014706 

PWALL 

0.065597 

PW/PO 

0.071217 

S 

1.919862 

X 

1.345333 

R 

1.054093 

NSH 

0.715520 

NSHP 

0.402819 

XSH 

1.251939 

RSH 

l.  763491 

NO  ITEP 
42 

USH 

VSH 

TSH 

RSH 

PSH 

0. 843984 

0.238788  0.129191 

5.373517 

0.198377 

USP 

VSP 

TSP 

RSP 

PSP 

0.106750 

-0.090286  -0.068657 

-0.331244 

-0.117652 

CF 

HEAT 

STAN 

COF 

CDP 

COTOT 

PWALL 

PW/PO 

033008  0. 

008852 

0.018634 

0.135843  0.857802  0.993646 

0.060939 

0.066214 

S 

X 

R 

NSH 

NSHP 

XSH 

RSH 

NO  ITER 

056955  l. 

481293 

1.071992 

0.769285  0.381310  1.380881 

1.834696 

42 

USH 

VSH 

TSH 

RSH 

PSH 

0.857710 

0.226726  0.120330 

5.327713 

0.183194 

USP 

VSP 

TSP 

RSP 

PSP 

0.091665 

-0.084696  -0.059486 

-0.330612 

-0. 101929 
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CF 

HEAT 

STAN 

CCF 

COP 

COTOT 

PWALL 

PW/PO 

0.030649 

0.008257 

0.017353 

0.139349 

0.833314 

0.972663 

0.057375 

0.062290 

s 

X 

R 

NSH 

NSHP 

XSH 

RSH 

NO  ITER 

2.  194128 

1.617253 

1.089892 

0.819539 

0.351615 

1.510282 

1.902419 

43 

USH 

VSH 

TSH 

RSH 

PSH 

0.869244 

0.215709  0.112817  5.283300  0 

.170320 

USP 

VSP 

TSP 

RSP 

PSP 

0.077149 

-0.077060  -0.050508  -0.319154  -0 

.086540 

CF 

HEAT 

STAN 

COF 

COP 

COTOT 

PWALL 

PK/PO 

0.028774 

0.007799 

0.016391 

0.142159 

0.809813 

0.951972 

0.054544 

0.059216 

s 

X 

R 

NSH 

NSHP 

XSH 

RSH 

NO  ITER 

2.331261 

1.753213 

1.107791 

0.864317 

0.301446 

1.640397 

1.964714 

45 

USH 

VSH 

TSH 

RSH 

DSH 

0.878704 

0.206006 

0.106604 

5.241906 

0.159675 

USP 

VSP 

TSP 

RSP 

PSP 

0.059373 

-0.063303 

-0.039178 

-0.277096 

-0.067122 

CF 

HEAT 

STAN 

CDF 

COP 

COTOT 

PWALL 

PW/PO 

0.027154 

0.007448 

0.015654 

0. 144409 

0.787287 

0.931696 

0.052554 

0.057056 

S 

X 

R 

NSH 

NSHP 

XSH 

RSH 

NO  ITER 

2. 468394 

1.889172 

1. 125691 

0.904733 

0.288002 

1.771081 

2.022684 

42 

USH 

VSH 

TSH 

RSH 

PSH 

0.886844 

0.197100 

0.101215 

5.201760 

0.150442 

USP 

VSP 

TSP 

RSP 

PSP 

0.050353 

-0.056688 

-0.033461 

-0.262395 

-0.057324 

CF 

HEAT 

STAN 

COF 

CDP 

COTOT 

PWALL 

PW/PO 

0.025854  0. 

007137 

0.014999 

0. 146203 

0.765712 

0.911916 

0.050898 

0.055258 

$ 

X 

R 

NSH 

NSHP 

XSH 

RSH 

NO  ITER 

2.  605527  2. 

025132 

1. 143590 

0.950673 

0.382001 

1.901045 

2.086130 

48 

USH 

VSH 

TSH 

RSH 

PSH 

0. 895073 

0.187512  0.095721  5.155987  0 

.141031 

USP 

VSP 

TSP 

RSP 

PSP 

0.067366 

-0.080597  -0.045135  -0.395463  -0 

.077315 

CF 

HEAT 

STAN 

CCF 

CCP 

COTOT 

PWALL 

PK/PO 

0.025167  0. 

006817 

0.014328 

0. 147680 

0.745016 

0. 892696 

0.048489 

0.052643 
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SYMBOLS 

c*  viscosity  law  constant,  c*  = 198. 6°R 

* * *2» 

Cf  skin  friction  coefficient,  2t  /(p^  u^  ) 

C*  specific  heat  of  constant  pressure 

^ * * 2 
H nondimensional  total  enthalpy,  H /u^ 

k thermal  conductivity 

M free  stream  Mach  number 

00 

n coordinate  measured  normal  to  the  body,  nondimen- 

sionalized  by  the  body  nose  radius 

n shock  stand  off  distance  normal  to  the  body  surface 

s 

* * 2*. 

p nondimensional  pressure,  p /( p ^ ub  ) 

* * * 3 

q nondimensional  heat  transfer,  q /(p^  ua  ) 

r nondimensional  axisymmetric  radius 

R defined  as  yB  + ng  costfr 

s nondimensional  surface  distance  coordinate 

st  Stanton  number,  st  = qw/(HQ-Hw) 

T nondimensional  temperature,  T = T /(.u^  /C  ) 

T*  free  stream  temperature 

00 

u nondimensional  velocity  component  tangent  to  .the 

body  surface,  u*/u*  . 

•fa  s 

u free  stream  velocity 

CO 

u nondimensional  component  of  velocity  aft  and 

tangent  to  the  shock  interface 

v nondimensional  velocity  component  normal  to  the 

body  surface,  v*/u* 

v nondimensional  component  of  velocity  aft  and 

normal  to  shock  interface 
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x„  axial  distance  for  body  surface  measured  from 
stagnation  point 

xg  defined  as  xg  - ng  sin<j> 

y normal  distance  for  body  surface  measured 

from  axis 

a shock  angle,  see  Figure  1 

0 angle  defined  in  Figure  1 , 

y ratio  of  specific  heats 

* *2  * * * * i/2 

e perturbation  parameter,  e = [y  (um  /C^J/p^u^a  ] ' 

k nondimensional  surface  curvature 

y nondimensional  coefficient  of  viscosity, 

* * *2  * 
y = y /y  (u^  /Cp) 

* * 

p nondimensional  density,  p = p /pa 

* 

p^  free  stream  density 

ic  'fc  ic  o 

r nondimensional  shear  stress,  t / ( p u ) 

00  00 

■ body  angle  defined  in  Figure  1 
a Prandtl  number,  a = yCp/K 

Subscripts 

1 wall  value 

0 stagnation  conditions 

s used  for  longitudinal  derivatives 

sh  conditions  immediately  behind  the  shock  wave 

» free  stream  conditions 
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Superscripts 


physical  quantities  normalized  by  their  shock  values 

* dimensional  quantities,  also  used  for  first  sweep 

of  ADI  numerical  scheme 

J 0 for  plane  flow  and  1 for  axisymmetric  flow 

n+1  represents  second  sweep  of  ADI  numerical  scheme 
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